From 3c043af4dbd8119fc8966c220352d4d9976e4317 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Fri, 11 Sep 2026 14:42:55 +0200 Subject: [PATCH] musrview: rewrite on top of the ROOT7 web graphics musrview is now a client/server application following the ROOT7 architecture (S. Linev et al., J. Phys. Conf. Ser. 1525 (2020) 012092): the C++ side holds the data, an ROOT::RWebWindow ships it as JSON to a JSROOT page which renders it. This addresses three long standing problems. The fit parameter, theory and run info blocks were TPaveText/TLegend objects which scaled their content down until it became unreadable as soon as a fit had many parameters. They are HTML panels now: fixed font size, own scrollbar, and each of them can be shown and hidden individually. The run info keeps listing every run with the very symbol it is plotted with. The command keys ('d', 'f', 'a', ...) were executed even while a text or latex label was being typed into the canvas. The page checks whether a text input has the focus, so typing a 'd' inserts a 'd'. The plot views were encoded in a set of entangled variables plus a static one which leaked state between canvases. PMusrViewMode describes them as four orthogonal properties - domain, Fourier component, difference including whether it is formed before or after the transform, and average - so that combinations such as 'Fourier-difference-average' are well defined. PMusrViewModel::GetCurves() is the single place which maps a view mode onto histograms; plotting and ascii export both go through it and can no longer disagree, which removed some 700 lines of duplicated view logic from the export. New classes, namespace musrview: PMusrViewMode the decoupled view state, free of ROOT and GUI, unit tested PMusrViewModel the GUI-free data layer, physics taken over unchanged PMusrWebCanvas the web front end and the graphics export Further notes: - The interactive path builds no TCanvas at all, the spectra travel as a TList of drawables. A TCanvas creates and maps its window already in the constructor, which would put an empty ROOT canvas next to the musrview window. - The graphics export keeps the classic pad layout and its 800x600 size and runs with ROOT in batch mode, so batch operation needs no browser. - The window opens at 1200x900, since toolbar and panels take their share. - Non-muSR axis labels are taken from the data file again, and ROOT/TLatex markup is rendered in the HTML panels as well. - musrview needs a ROOT built with -Droot7=ON -Dwebgui=ON, plus -Dqt6web=ON for the embedded window (--embedded, the default). --browser uses the standard browser of the system, --browser=server only prints the URL. - PMusrCanvas is no longer part of the build, the files are kept for reference. Tests: 87/87 pass, among them a new unit test of the view modes and a new musrview test covering the non-muSR fit type. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01H6YiUNBTMDX9uA3Ts1w2nW --- CLAUDE.md | 32 +- CMakeLists.txt | 10 + ChangeLog | 50 + src/CMakeLists.txt | 10 +- src/classes/CMakeLists.txt | 30 +- src/classes/PMusrViewMode.cpp | 486 ++ src/classes/PMusrViewModel.cpp | 3914 +++++++++++++++++ src/classes/PMusrWebCanvas.cpp | 1117 +++++ src/include/PMusrViewMode.h | 199 + src/include/PMusrViewModel.h | 391 ++ src/include/PMusrWebCanvas.h | 185 + src/musrview.cpp | 346 +- src/musrview_web/musrview.html | 705 +++ tests/CMakeLists.txt | 1 + tests/musrview_check/CMakeLists.txt | 29 + .../test-libGapIntegrals-ASCII_0.png | Bin 0 -> 17383 bytes tests/viewmode_check/CMakeLists.txt | 15 + tests/viewmode_check/viewmode_check.cpp | 88 + 18 files changed, 7460 insertions(+), 148 deletions(-) create mode 100644 src/classes/PMusrViewMode.cpp create mode 100644 src/classes/PMusrViewModel.cpp create mode 100644 src/classes/PMusrWebCanvas.cpp create mode 100644 src/include/PMusrViewMode.h create mode 100644 src/include/PMusrViewModel.h create mode 100644 src/include/PMusrWebCanvas.h create mode 100644 src/musrview_web/musrview.html create mode 100644 tests/musrview_check/ref/musrview-nonMusr-GapIntegrals/test-libGapIntegrals-ASCII_0.png create mode 100644 tests/viewmode_check/CMakeLists.txt create mode 100644 tests/viewmode_check/viewmode_check.cpp diff --git a/CLAUDE.md b/CLAUDE.md index 6839d3ea..5cb6a125 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -92,8 +92,36 @@ The central data model and control flow (all in `src/classes/` with headers in ` in the msr file). - **`PFitterFcn`** wraps the χ²/max-log-likelihood objective for ROOT's Minuit2; **`PFitter`** drives the minimization and post-fit error analysis. -- **`PMusrCanvas`** / **`PFourier`** / **`PFourierCanvas`** / **`PPrepFourier`** handle plotting - and Fourier-transform display, used by `musrview`/`musrFT`. +- **`PMusrViewMode`** / **`PMusrViewModel`** / **`PMusrWebCanvas`** (namespace `musrview`) are + what `musrview` is built from since the ROOT7 rewrite: + - `PMusrViewMode` describes what is shown as four *orthogonal* properties: domain + (time/Fourier), Fourier component, difference (and whether it is formed before or after the + transform, d-f vs f-d) and average. No GUI/ROOT dependency; unit tested in + `tests/viewmode_check/`. + - `PMusrViewModel` is the GUI-free data layer: it builds the data/theory histograms, their + Fourier transforms, the differences and the average, and answers `GetCurves(mode)` with the + histograms that mode makes visible. Plotting and ascii export both go through that one + selector, so they cannot disagree. + - `PMusrWebCanvas` is the front end. `CollectDrawables()` returns the objects of the current + view with their ROOT draw options; from that, the interactive path serializes a plain `TList` + with `TBufferJSON` and pushes it through an `ROOT::RWebWindow` to a JSROOT page + (`src/musrview_web/musrview.html`), which draws them in order and shows + title/parameters/theory/run-info as scrollable, individually collapsible HTML panels. + Deliberately no `TCanvas` is built for that path: a TCanvas maps its window in the + constructor, which would put an empty ROOT window next to the musrview window (toggling + `gROOT->SetBatch()` at run time crashes in `TCanvas::Build()`, and ROOT web canvases open + their own window). The graphics export (`--png`, ...) uses the same collector but draws into + a real `TCanvas` with the classic pad layout, running with ROOT in batch mode, so batch runs + need no browser. + - The web assets are installed to `/share/musrfit/musrview`; set `MUSRFIT_WEB_DIR` to run + from the source tree, and `MUSRVIEW_DEBUG=1` to trace the websocket traffic. + - musrview needs a ROOT built with `-Droot7=ON -Dwebgui=ON`, plus `-Dqt6web=ON` for the embedded + window (`--embedded`, the default). `--browser` uses the system browser, `--browser=server` + only prints the URL. +- **`PFourier`** / **`PFourierCanvas`** / **`PPrepFourier`** handle the Fourier transform and its + (still TCanvas based) display, used by `musrview`/`musrFT`. +- **`PMusrCanvas`** is the former TCanvas based musrview front end. It is no longer part of the + build; the files are kept for reference only. - **`PMsr2Data`** implements `msr2data`'s logic for generating follow-up/global msr files from a template across a run list. diff --git a/CMakeLists.txt b/CMakeLists.txt index 56fbee94..5371a90d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -101,6 +101,16 @@ set(MUSRFIT_CXX_FLAGS_NO_ROOT "${CMAKE_CXX_FLAGS}") find_package(ROOT 6.18 REQUIRED COMPONENTS Gui MathMore Minuit2 XMLParser) +#--- musrview is based on the ROOT7 web graphics, which needs a ROOT built +#--- with -Dwebgui=ON (and -Droot7=ON). Without it there is no RWebWindow and +#--- musrview cannot be built. ------------------------------------------------ +if (NOT TARGET ROOT::ROOTWebDisplay) + message(FATAL_ERROR + "musrview needs the ROOT web display (RWebWindow), which this ROOT " + "installation does not provide. Please rebuild ROOT with " + "-Droot7=ON -Dwebgui=ON (and -Dqt6web=ON for the embedded window).") +endif () + #--- Deliberately NOT 'include(${ROOT_USE_FILE})' here. That file #--- - overwrites CMAKE_CXX_STANDARD with ROOT's value, silently discarding #--- the setting made at the top of this file, diff --git a/ChangeLog b/ChangeLog index 2c38e634..9352fb47 100644 --- a/ChangeLog +++ b/ChangeLog @@ -12,6 +12,56 @@ or https://bitbucket.org/muonspin/musrfit +Unreleased +========== + +musrview has been rewritten on top of the ROOT7 web graphics (RWebWindow, +JSROOT), following the architecture described in S. Linev et al., +J. Phys. Conf. Ser. 1525 (2020) 012092. The C++ side is now a pure data server, +the presentation happens in a web page. This addresses three long standing +shortcomings: + +- The fit parameter, theory and run info blocks used to be TPaveText/TLegend + objects which scaled their content down until it became unreadable as soon as + a fit had many parameters. They are now HTML panels with a fixed font size and + their own scrollbar, and each of them can be shown and hidden individually. + The run info keeps listing every run with the very symbol it is plotted with, + colour and marker shape included, as the TLegend did. +- The command keys ('d', 'f', 'a', ...) were executed even while a text or latex + label was being typed into the canvas. The page now checks whether a text + input has the focus, so typing a 'd' inserts a 'd'. +- The plot views (time/Fourier, difference before or after the Fourier + transform, average) were encoded in a set of entangled variables plus a static + one which leaked state between canvases. They are now four orthogonal + properties of the new class PMusrViewMode, so that combinations such as + 'Fourier-difference-average' are well defined. The ascii export no longer + repeats the view logic but asks the model what is currently shown, which means + export and display can no longer disagree. + +The musrview window opens larger than the ROOT6 canvas used to (1200x900 +instead of 800x600), since the toolbar and the panels take their share of it. +The graphics export is unaffected and keeps producing 800x600 images. + +musrview can either open its window in the standard browser of the system +(--browser) or in an embedded Qt6 WebEngine window (--embedded, the default), +in which case it appears as the stand-alone application it always was. Use +--browser=server when musrfit runs on a remote machine; it only prints the URL. +This needs a ROOT built with -Droot7=ON -Dwebgui=ON, and -Dqt6web=ON for the +embedded window. + +The graphics export (--png, --pdf, ...) deliberately does not go through a +browser: it keeps producing the classic canvas layout via TCanvas, so batch +operation does not depend on a browser installation. + +In the interactive case no TCanvas is built at all. The spectra travel to +JSROOT as a plain TList of drawables, each with the draw option it has to be +drawn with. A TCanvas creates and maps its window already in the constructor, +so an empty ROOT canvas window would appear next to the musrview window, and +the ways around that are worse: switching gROOT to batch at run time leaves +ROOT's graphics state inconsistent and crashes in TCanvas::Build(), and ROOT +web canvases open a window of their own. Sending the drawables is also smaller +on the wire and is what the ROOT7 client/server split suggests anyway. + Release of V1.12.1, 2026/08/25 ============================== diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c097a7ec..58c7058a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -135,7 +135,15 @@ target_include_directories(musrview $ $ ) -target_link_libraries(musrview FFTW3::FFTW3 ${ROOT_LIBRARIES} ${MUSRFIT_LIBS} Boost::headers) +target_link_libraries(musrview FFTW3::FFTW3 ${ROOT_LIBRARIES} ROOT::ROOTWebDisplay ${MUSRFIT_LIBS} Boost::headers) + +#--- musrview is web based: the page it serves has to be installed as well ---- +#--- MUSRVIEW_WEB_DIR (see src/classes/CMakeLists.txt) points here. During +#--- development the environment variable MUSRFIT_WEB_DIR overrides it. ------- +install( + FILES ${CMAKE_SOURCE_DIR}/src/musrview_web/musrview.html + DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/musrfit/musrview +) add_executable(write_musrRoot_runHeader write_musrRoot_runHeader.cpp) target_compile_options(write_musrRoot_runHeader BEFORE PRIVATE "-DHAVE_CONFIG_H" "${HAVE_GIT_REV_H}") diff --git a/src/classes/CMakeLists.txt b/src/classes/CMakeLists.txt index cb992e62..4bfa6fff 100644 --- a/src/classes/CMakeLists.txt +++ b/src/classes/CMakeLists.txt @@ -32,15 +32,6 @@ root_generate_dictionary( -inlineInputHeader MODULE PMsgBox ) -root_generate_dictionary( - PMusrCanvasDict - PMusrCanvas.h - LINKDEF ${MUSRFIT_INC}/PMusrCanvasLinkDef.h - OPTIONS - -isystem ${Boost_INCLUDE_DIRS} -isystem ${FFTW3_INCLUDE} -I${MUSRFIT_INC} - -inlineInputHeader - MODULE PMusrCanvas -) root_generate_dictionary( PMusrT0Dict PMusrT0.h @@ -103,8 +94,9 @@ add_library(PMusr SHARED PFourier.cpp PFourierCanvas.cpp PFourierCanvasDict.cxx - PMusrCanvas.cpp - PMusrCanvasDict.cxx + PMusrViewMode.cpp + PMusrViewModel.cpp + PMusrWebCanvas.cpp PFunction.cpp PFunctionHandler.cpp PMsgBox.cpp @@ -112,8 +104,6 @@ add_library(PMusr SHARED PMsr2Data.cpp PMsrHandler.cpp PStringUtils.cpp - PMusrCanvas.cpp - PMusrCanvasDict.cxx PMusr.cpp PMusrT0.cpp PMusrT0Dict.cxx @@ -160,6 +150,12 @@ add_library(PRgeHandler SHARED #--- set target properties, e.g. version -------------------------------------- target_compile_options(PMusr BEFORE PRIVATE "-DHAVE_CONFIG_H") +#--- where musrview finds its web assets (html page) at runtime --------------- +#--- MUSRVIEW_WEB_DIR is the installed location, MUSRVIEW_WEB_DIR_DEV the +#--- source tree so that musrview also runs straight from the build directory. +target_compile_definitions(PMusr PRIVATE + MUSRVIEW_WEB_DIR="${CMAKE_INSTALL_PREFIX}/${CMAKE_INSTALL_DATAROOTDIR}/musrfit/musrview" + MUSRVIEW_WEB_DIR_DEV="${CMAKE_SOURCE_DIR}/src/musrview_web") set_target_properties(PMusr PROPERTIES VERSION ${MUSR_VERSION} @@ -197,6 +193,8 @@ endif(ROOT_GRTEQ_24) set(DependOnLibs ${ROOT_LIBRARIES}) #set(DependOnLibs ${DependOnLibs} ${FFTW3_LIBRARY}) set(DependOnLibs ${DependOnLibs} FFTW3::FFTW3) +#--- musrview (PMusrWebCanvas) is built on ROOT7 RWebWindow ------------------- +set(DependOnLibs ${DependOnLibs} ROOT::ROOTWebDisplay) set(DependOnLibs ${DependOnLibs} PUserFcnBase) set(DependOnLibs ${DependOnLibs} mud) set(DependOnLibs ${DependOnLibs} TMusrRunHeader) @@ -230,8 +228,6 @@ install( ${CMAKE_CURRENT_BINARY_DIR}/libPFourierCanvas.rootmap ${CMAKE_CURRENT_BINARY_DIR}/libPMsgBox_rdict.pcm ${CMAKE_CURRENT_BINARY_DIR}/libPMsgBox.rootmap - ${CMAKE_CURRENT_BINARY_DIR}/libPMusrCanvas_rdict.pcm - ${CMAKE_CURRENT_BINARY_DIR}/libPMusrCanvas.rootmap ${CMAKE_CURRENT_BINARY_DIR}/libPMusrT0_rdict.pcm ${CMAKE_CURRENT_BINARY_DIR}/libPMusrT0.rootmap ${CMAKE_CURRENT_BINARY_DIR}/libPRgeHandler_rdict.pcm @@ -257,8 +253,10 @@ install( ${MUSRFIT_INC}/PMsgBox.h ${MUSRFIT_INC}/PMsr2Data.h ${MUSRFIT_INC}/PMsrHandler.h - ${MUSRFIT_INC}/PMusrCanvas.h ${MUSRFIT_INC}/PMusr.h + ${MUSRFIT_INC}/PMusrViewMode.h + ${MUSRFIT_INC}/PMusrViewModel.h + ${MUSRFIT_INC}/PMusrWebCanvas.h ${MUSRFIT_INC}/PMusrT0.h ${MUSRFIT_INC}/PPrepFourier.h ${MUSRFIT_INC}/PRgeHandler.h diff --git a/src/classes/PMusrViewMode.cpp b/src/classes/PMusrViewMode.cpp new file mode 100644 index 00000000..44389ac6 --- /dev/null +++ b/src/classes/PMusrViewMode.cpp @@ -0,0 +1,486 @@ +/*************************************************************************** + + PMusrViewMode.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2026 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. * + ***************************************************************************/ + +#include "PMusrViewMode.h" +#include "PMusr.h" + +namespace musrview { + +//-------------------------------------------------------------------------- +// SetFourierAllowed (public) +//-------------------------------------------------------------------------- +/** + *

Allows/forbids the Fourier domain. If the Fourier domain is forbidden while + * it is the current one, the view falls back to the time domain. + * + * \param allowed true if the data can be Fourier transformed + */ +void PMusrViewMode::SetFourierAllowed(bool allowed) +{ + fFourierAllowed = allowed; + if (!fFourierAllowed) + fDomain = EDomain::kTime; +} + +//-------------------------------------------------------------------------- +// SetAverageAllowed (public) +//-------------------------------------------------------------------------- +/** + *

Allows/forbids averaging. If averaging is forbidden while it is active, it + * is switched off. + * + * \param allowed true if the data sets can be averaged + */ +void PMusrViewMode::SetAverageAllowed(bool allowed) +{ + fAverageAllowed = allowed; + if (!fAverageAllowed) + fAverage = false; +} + +//-------------------------------------------------------------------------- +// ToggleFourier (public) +//-------------------------------------------------------------------------- +/** + *

Switches between time and Fourier domain, the 'f' command key. + * + *

If the difference is active while switching into the Fourier domain, the + * difference already exists as a time spectrum, hence it is this difference + * which gets transformed (d-f). + * + * \return true if the state has changed + */ +bool PMusrViewMode::ToggleFourier() +{ + return SetDomain(IsFourier() ? EDomain::kTime : EDomain::kFourier); +} + +//-------------------------------------------------------------------------- +// ToggleDifference (public) +//-------------------------------------------------------------------------- +/** + *

Switches the difference data-theory on/off, the 'd' command key. + * + * \return true if the state has changed + */ +bool PMusrViewMode::ToggleDifference() +{ + return SetDifference(!fDifference); +} + +//-------------------------------------------------------------------------- +// ToggleAverage (public) +//-------------------------------------------------------------------------- +/** + *

Switches averaging on/off, the 'a' command key. + * + * \return true if the state has changed + */ +bool PMusrViewMode::ToggleAverage() +{ + return SetAverage(!fAverage); +} + +//-------------------------------------------------------------------------- +// SetDomain (public) +//-------------------------------------------------------------------------- +/** + *

Selects the representation domain. + * + *

Entering the Fourier domain while a difference is shown means that the + * time domain difference gets transformed, i.e. d-f. + * + * \param domain the requested domain + * + * \return true if the state has changed + */ +bool PMusrViewMode::SetDomain(EDomain domain) +{ + if ((domain == EDomain::kFourier) && !fFourierAllowed) + return false; + if (domain == fDomain) + return false; + + if (domain == EDomain::kFourier) + fDiffBeforeFourier = true; // the existing time domain difference is transformed + + fDomain = domain; + + return true; +} + +//-------------------------------------------------------------------------- +// SetDifference (public) +//-------------------------------------------------------------------------- +/** + *

Switches the difference data-theory on/off. + * + *

Switching the difference on while Fourier spectra are shown means that the + * Fourier spectra of data and theory are subtracted, i.e. f-d. + * + * \param on true if the difference shall be shown + * + * \return true if the state has changed + */ +bool PMusrViewMode::SetDifference(bool on) +{ + if (on == fDifference) + return false; + + if (on) + fDiffBeforeFourier = IsTimeDomain(); + + fDifference = on; + + return true; +} + +//-------------------------------------------------------------------------- +// SetAverage (public) +//-------------------------------------------------------------------------- +/** + *

Switches averaging over all data sets on/off. + * + * \param on true if the average shall be shown + * + * \return true if the state has changed + */ +bool PMusrViewMode::SetAverage(bool on) +{ + if (on && !fAverageAllowed) + return false; + if (on == fAverage) + return false; + + fAverage = on; + + return true; +} + +//-------------------------------------------------------------------------- +// SetFourierComponent (public) +//-------------------------------------------------------------------------- +/** + *

Selects which part of the complex Fourier transform is shown. This implies + * the Fourier domain. + * + * \param comp the requested Fourier component + * + * \return true if the state has changed + */ +bool PMusrViewMode::SetFourierComponent(EFourierComponent comp) +{ + if (!fFourierAllowed) + return false; + + bool changed = SetDomain(EDomain::kFourier); + + if (comp != fFourierComponent) { + fFourierComponent = comp; + changed = true; + } + + return changed; +} + +//-------------------------------------------------------------------------- +// Reset (public) +//-------------------------------------------------------------------------- +/** + *

Resets the view to the plain time domain representation. The capabilities + * (Fourier/average allowed) are kept. + */ +void PMusrViewMode::Reset() +{ + fDomain = EDomain::kTime; + fDifference = false; + fDiffBeforeFourier = true; + fAverage = false; +} + +//-------------------------------------------------------------------------- +// GetDiffStage (public) +//-------------------------------------------------------------------------- +/** + *

Where the difference is formed relative to the Fourier transform. + * + * \return EDiffStage::kNone if no difference is shown, otherwise the stage. In + * the time domain the placement is meaningless and the canonical value + * EDiffStage::kBeforeFourier is returned. + */ +EDiffStage PMusrViewMode::GetDiffStage() const +{ + if (!fDifference) + return EDiffStage::kNone; + + if (IsTimeDomain()) + return EDiffStage::kBeforeFourier; + + return fDiffBeforeFourier ? EDiffStage::kBeforeFourier : EDiffStage::kAfterFourier; +} + +//-------------------------------------------------------------------------- +// IsPhaseAdjustable (public) +//-------------------------------------------------------------------------- +/** + *

A Fourier phase can only be applied to a complex spectrum. The power + * spectrum is invariant under a phase rotation, hence it is excluded. + * + * \return true if the '+'/'-' phase commands have an effect + */ +bool PMusrViewMode::IsPhaseAdjustable() const +{ + return IsFourier() && (fFourierComponent != EFourierComponent::kPower); +} + +//-------------------------------------------------------------------------- +// GetFourierComponentName (static public) +//-------------------------------------------------------------------------- +/** + *

Machine readable name of a Fourier component. + * + * \param comp Fourier component + * + * \return the name, e.g. "real" + */ +std::string PMusrViewMode::GetFourierComponentName(EFourierComponent comp) +{ + switch (comp) { + case EFourierComponent::kReal: return "real"; + case EFourierComponent::kImag: return "imag"; + case EFourierComponent::kRealAndImag: return "real-imag"; + case EFourierComponent::kPower: return "power"; + case EFourierComponent::kPhase: return "phase"; + case EFourierComponent::kPhaseOptReal: return "phase-opt-real"; + } + return "power"; +} + +//-------------------------------------------------------------------------- +// GetFourierComponentLabel (static public) +//-------------------------------------------------------------------------- +/** + *

Human readable name of a Fourier component. + * + * \param comp Fourier component + * + * \return the label, e.g. "real part" + */ +std::string PMusrViewMode::GetFourierComponentLabel(EFourierComponent comp) +{ + switch (comp) { + case EFourierComponent::kReal: return "real part"; + case EFourierComponent::kImag: return "imaginary part"; + case EFourierComponent::kRealAndImag: return "real+imaginary part"; + case EFourierComponent::kPower: return "power"; + case EFourierComponent::kPhase: return "phase"; + case EFourierComponent::kPhaseOptReal: return "phase optimized real part"; + } + return "power"; +} + +//-------------------------------------------------------------------------- +// FourierComponentFromMsrTag (static public) +//-------------------------------------------------------------------------- +/** + *

Maps the FOURIER_PLOT_* tag found in the msr-file FOURIER block onto the + * corresponding Fourier component. + * + * \param plotTag FOURIER_PLOT_* tag + * + * \return the Fourier component; the power spectrum is used as fallback + */ +EFourierComponent PMusrViewMode::FourierComponentFromMsrTag(int plotTag) +{ + switch (plotTag) { + case FOURIER_PLOT_REAL: return EFourierComponent::kReal; + case FOURIER_PLOT_IMAG: return EFourierComponent::kImag; + case FOURIER_PLOT_REAL_AND_IMAG: return EFourierComponent::kRealAndImag; + case FOURIER_PLOT_POWER: return EFourierComponent::kPower; + case FOURIER_PLOT_PHASE: return EFourierComponent::kPhase; + case FOURIER_PLOT_PHASE_OPT_REAL: return EFourierComponent::kPhaseOptReal; + default: return EFourierComponent::kPower; + } +} + +//-------------------------------------------------------------------------- +// MsrTagFromFourierComponent (static public) +//-------------------------------------------------------------------------- +/** + *

Inverse of FourierComponentFromMsrTag(). + * + * \param comp Fourier component + * + * \return the corresponding FOURIER_PLOT_* tag + */ +int PMusrViewMode::MsrTagFromFourierComponent(EFourierComponent comp) +{ + switch (comp) { + case EFourierComponent::kReal: return FOURIER_PLOT_REAL; + case EFourierComponent::kImag: return FOURIER_PLOT_IMAG; + case EFourierComponent::kRealAndImag: return FOURIER_PLOT_REAL_AND_IMAG; + case EFourierComponent::kPower: return FOURIER_PLOT_POWER; + case EFourierComponent::kPhase: return FOURIER_PLOT_PHASE; + case EFourierComponent::kPhaseOptReal: return FOURIER_PLOT_PHASE_OPT_REAL; + } + return FOURIER_PLOT_POWER; +} + +//-------------------------------------------------------------------------- +// GetKey (public) +//-------------------------------------------------------------------------- +/** + *

Stable machine readable representation of the state, e.g. + * "fourier-real+diff-fd+avg". + * + * \return the key + */ +std::string PMusrViewMode::GetKey() const +{ + std::string key; + + if (IsFourier()) + key = "fourier-" + GetFourierComponentName(fFourierComponent); + else + key = "time"; + + switch (GetDiffStage()) { + case EDiffStage::kNone: + break; + case EDiffStage::kBeforeFourier: + key += IsFourier() ? "+diff-df" : "+diff"; + break; + case EDiffStage::kAfterFourier: + key += "+diff-fd"; + break; + } + + if (fAverage) + key += "+avg"; + + return key; +} + +//-------------------------------------------------------------------------- +// GetLabel (public) +//-------------------------------------------------------------------------- +/** + *

Human readable description of the state, used as plot title decoration. + * + * \return the label + */ +std::string PMusrViewMode::GetLabel() const +{ + std::string label; + + if (IsFourier()) + label = "Fourier " + GetFourierComponentLabel(fFourierComponent); + else + label = "Data"; + + switch (GetDiffStage()) { + case EDiffStage::kNone: + break; + case EDiffStage::kBeforeFourier: + label += IsFourier() ? ", difference (d-f)" : ", difference"; + break; + case EDiffStage::kAfterFourier: + label += ", difference (f-d)"; + break; + } + + if (fAverage) + label += ", averaged"; + + return label; +} + +//-------------------------------------------------------------------------- +// ToJson (public) +//-------------------------------------------------------------------------- +/** + *

JSON representation of the state which is sent to the web front end so + * that it can render the menu/toolbar without duplicating any of the logic. + * + * \return the JSON object as string + */ +std::string PMusrViewMode::ToJson() const +{ + std::string diffStage; + switch (GetDiffStage()) { + case EDiffStage::kNone: diffStage = "none"; break; + case EDiffStage::kBeforeFourier: diffStage = "df"; break; + case EDiffStage::kAfterFourier: diffStage = "fd"; break; + } + + std::string json = "{"; + json += "\"key\":\"" + GetKey() + "\","; + json += "\"label\":\"" + GetLabel() + "\","; + json += "\"domain\":\"" + std::string(IsFourier() ? "fourier" : "time") + "\","; + json += "\"fourierComponent\":\"" + GetFourierComponentName(fFourierComponent) + "\","; + json += "\"difference\":" + std::string(fDifference ? "true" : "false") + ","; + json += "\"diffStage\":\"" + diffStage + "\","; + json += "\"average\":" + std::string(fAverage ? "true" : "false") + ","; + json += "\"fourierAllowed\":" + std::string(fFourierAllowed ? "true" : "false") + ","; + json += "\"averageAllowed\":" + std::string(fAverageAllowed ? "true" : "false") + ","; + json += "\"phaseAdjustable\":" + std::string(IsPhaseAdjustable() ? "true" : "false"); + json += "}"; + + return json; +} + +//-------------------------------------------------------------------------- +// operator== (public) +//-------------------------------------------------------------------------- +/** + *

Two modes are equal if they select the same representation. The + * capabilities are properties of the data, not of the view, and are therefore + * not compared. + * + * \param rhs mode to compare with + * + * \return true if both modes show the same thing + */ +bool PMusrViewMode::operator==(const PMusrViewMode &rhs) const +{ + if (fDomain != rhs.fDomain) + return false; + if (fAverage != rhs.fAverage) + return false; + if (fDifference != rhs.fDifference) + return false; + if (GetDiffStage() != rhs.GetDiffStage()) + return false; + if (IsFourier() && (fFourierComponent != rhs.fFourierComponent)) + return false; + + return true; +} + +} // namespace musrview diff --git a/src/classes/PMusrViewModel.cpp b/src/classes/PMusrViewModel.cpp new file mode 100644 index 00000000..c95fdc5c --- /dev/null +++ b/src/classes/PMusrViewModel.cpp @@ -0,0 +1,3914 @@ +/*************************************************************************** + + PMusrViewModel.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2026 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 +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "PMusrViewModel.h" +#include "PFourier.h" +#include "PStringUtils.h" + +namespace musrview { + +//-------------------------------------------------------------------------- +// PMusrViewPlotRange +//-------------------------------------------------------------------------- +/** + *

Sets the x-range and marks it as present. The values are swapped if given + * in the wrong order. + * + * \param xmin minimum x value + * \param xmax maximum x value + */ +void PMusrViewPlotRange::SetXRange(Double_t xmin, Double_t xmax) +{ + if (xmin > xmax) { + std::cerr << std::endl << ">> PMusrViewPlotRange::SetXRange(): **WARNING** xmin > xmax, will swap them." << std::endl; + fXmin = xmax; + fXmax = xmin; + } else { + fXmin = xmin; + fXmax = xmax; + } + fXRangePresent = true; +} + +//-------------------------------------------------------------------------- +/** + *

Sets the y-range and marks it as present. The values are swapped if given + * in the wrong order. + * + * \param ymin minimum y value + * \param ymax maximum y value + */ +void PMusrViewPlotRange::SetYRange(Double_t ymin, Double_t ymax) +{ + if (ymin > ymax) { + std::cerr << std::endl << ">> PMusrViewPlotRange::SetYRange(): **WARNING** ymin > ymax, will swap them." << std::endl; + fYmin = ymax; + fYmax = ymin; + } else { + fYmin = ymin; + fYmax = ymax; + } + fYRangePresent = true; +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Creates the data layer for one PLOT block of the msr-file. + * + * \param plotNumber number of the PLOT block + * \param fourierDefault Fourier settings from the startup file + * \param markerList list of markers from the startup file + * \param colorList list of colors from the startup file + * \param theoAsData if true, the theory is only calculated at the data points + */ +PMusrViewModel::PMusrViewModel(Int_t plotNumber, PMsrFourierStructure fourierDefault, + PIntVector markerList, PIntVector colorList, Bool_t theoAsData) : + fTheoAsData(theoAsData), fPlotNumber(plotNumber), fFourier(fourierDefault), + fMarkerList(std::move(markerList)), fColorList(std::move(colorList)) +{ + InitAverage(); +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Releases all histograms/graphs owned by the model. + */ +PMusrViewModel::~PMusrViewModel() +{ + for (auto &ds : fData) + CleanupDataSet(ds); + fData.clear(); + + for (auto &ds : fNonMusrData) + CleanupDataSet(ds); + fNonMusrData.clear(); + + CleanupAverage(); + + if (fHistoFrame) { + delete fHistoFrame; + fHistoFrame = nullptr; + } +} + +//-------------------------------------------------------------------------- +// BuildData (public) +//-------------------------------------------------------------------------- +/** + *

Builds all time domain data/theory histograms of this PLOT block, the + * frame histogram which defines the visible range, and the text blocks. + * + *

This is the former PMusrCanvas::UpdateDataTheoryPad() with all plotting + * removed. + * + * \return true if everything could be built + */ +Bool_t PMusrViewModel::BuildData() +{ + if (fMsrHandler == nullptr || fRunList == nullptr) { + fLastError = "msr-handler and/or run list not set."; + return false; + } + + PMsrPlotStructure plotInfo = fMsrHandler->GetMsrPlotList()->at(fPlotNumber); + PMsrRunList runs = *fMsrHandler->GetMsrRunList(); + PMsrGlobalBlock *globalBlock = fMsrHandler->GetMsrGlobal(); + + Int_t fitType = globalBlock->GetFitType(); + + fPlotType = plotInfo.fPlotType; + + // check that the plot block is consistent with the run/global blocks + for (UInt_t i=0; i(plotInfo.fRuns[i]) > static_cast(runs.size())) { + std::ostringstream err; + err << "run plot number " << static_cast(plotInfo.fRuns[i]) + << " is larger than the number of runs " << runs.size(); + fLastError = err.str(); + std::cerr << std::endl << ">> PMusrViewModel::BuildData(): **ERROR** " << fLastError << std::endl; + return false; + } + UInt_t runNo = static_cast(plotInfo.fRuns[i])-1; + if (runs[runNo].GetFitType() != -1) + fitType = runs[runNo].GetFitType(); + if (fitType == -1) { + std::ostringstream err; + err << "plottype = " << fPlotType << ", fittype = " << runs[runNo].GetFitType() + << " (RUN block)/fittype = " << globalBlock->GetFitType() + << " (GLOBAL block). However, they have to correspond!"; + fLastError = err.str(); + std::cerr << std::endl << ">> PMusrViewModel::BuildData(): **ERROR** " << fLastError << std::endl; + return false; + } + } + + // collect the data of all runs of this plot block + for (UInt_t i=0; i(plotInfo.fRuns[i])-1; + if (runs[runNo].GetFitType() != -1) + fitType = runs[runNo].GetFitType(); + + const char *what = ""; + switch (fitType) { + case MSR_FITTYPE_SINGLE_HISTO: + data = fRunList->GetSingleHisto(runNo, PRunListCollection::kRunNo); + what = "a single histogram plot"; + break; + case MSR_FITTYPE_SINGLE_HISTO_RRF: + data = fRunList->GetSingleHistoRRF(runNo, PRunListCollection::kRunNo); + what = "a single histogram RRF plot"; + break; + case MSR_FITTYPE_ASYM: + data = fRunList->GetAsymmetry(runNo, PRunListCollection::kRunNo); + what = "an asymmetry plot"; + break; + case MSR_FITTYPE_BNMR: + data = fRunList->GetAsymmetryBNMR(runNo, PRunListCollection::kRunNo); + what = "a beta-NMR asymmetry plot"; + break; + case MSR_FITTYPE_ASYM_RRF: + data = fRunList->GetAsymmetryRRF(runNo, PRunListCollection::kRunNo); + what = "an asymmetry RRF plot"; + break; + case MSR_FITTYPE_MU_MINUS: + data = fRunList->GetMuMinus(runNo, PRunListCollection::kRunNo); + what = "a mu minus single histogram plot"; + break; + case MSR_FITTYPE_NON_MUSR: + data = fRunList->GetNonMusr(runNo, PRunListCollection::kRunNo); + what = "a non-muSR plot"; + break; + default: + fLastError = "wrong plottype tag?!"; + std::cerr << std::endl << ">> PMusrViewModel::BuildData(): **ERROR** " << fLastError << std::endl; + return false; + } + + if (data == nullptr) { + std::ostringstream err; + err << "couldn't obtain run no " << runNo << " for " << what; + fLastError = err.str(); + std::cerr << std::endl << ">> PMusrViewModel::BuildData(): **ERROR** " << fLastError << std::endl; + return false; + } + + if (fitType == MSR_FITTYPE_NON_MUSR) + HandleNonMusrDataSet(i, runNo, data); + else + HandleDataSet(i, runNo, data); + } + + fValid = true; + + // the frame histogram defines the visible time range and hence also the + // window of the Fourier transform. It used to be created as a side effect of + // TPad::DrawFrame(), here the model owns it. + if (!IsNonMusr() && !fData.empty()) { + Double_t yMin = 0.0, yMax = 0.0; + GetDefaultYRange(PMusrViewMode(), fXmin, fXmax, yMin, yMax); + + UInt_t noOfPoints = 1000; + for (const auto &ds : fData) { + if (ds.data && ds.data->GetNbinsX() > static_cast(noOfPoints)) + noOfPoints = ds.data->GetNbinsX(); + } + noOfPoints *= 2; // make sure that there are enough points + + if (fHistoFrame) + delete fHistoFrame; + fHistoFrame = new TH1F("hframe", "hframe", noOfPoints, fXmin, fXmax); + fHistoFrame->SetDirectory(nullptr); + fHistoFrame->SetStats(kFALSE); + fHistoFrame->SetMinimum(yMin); + fHistoFrame->SetMaximum(yMax); + } + + BuildTextBlocks(); + + return true; +} + +//-------------------------------------------------------------------------- +// SetVisibleTimeRange (public) +//-------------------------------------------------------------------------- +/** + *

Informs the model about the time range the user currently looks at. The + * Fourier transform is calculated over exactly this window, which is what makes + * 'zoom in, then press f' behave as expected. + * + * \param xmin lower bound of the visible time range + * \param xmax upper bound of the visible time range + */ +void PMusrViewModel::SetVisibleTimeRange(Double_t xmin, Double_t xmax) +{ + if (fHistoFrame == nullptr) + return; + if (xmin >= xmax) + return; + + fHistoFrame->GetXaxis()->SetRangeUser(xmin, xmax); + fVisibleRangeSet = true; + + // the Fourier transforms were calculated for a different window + CleanupFourier(); + CleanupFourierDifference(); + fFourierPrepared = false; + fDiffFourierPrepared = false; + fFourierDiffPrepared = false; + fAveragePrepared = false; +} + +//-------------------------------------------------------------------------- +// CreateDefaultMode (public) +//-------------------------------------------------------------------------- +/** + *

Creates the initial view mode: it carries the capabilities of the loaded + * data and honours the FOURIER block of the msr-file as well as the command + * line switches. + * + * \param startWithFourier start directly with the Fourier representation + * \param startWithAvg start directly with the averaged representation + * + * \return the initial view mode + */ +PMusrViewMode PMusrViewModel::CreateDefaultMode(Bool_t startWithFourier, Bool_t startWithAvg) const +{ + PMusrViewMode mode; + + mode.SetFourierAllowed(!IsNonMusr()); + mode.SetAverageAllowed(GetNoOfDataSets() > 1); + + mode.SetFourierComponent(PMusrViewMode::FourierComponentFromMsrTag(fFourier.fPlotTag)); + mode.SetDomain(EDomain::kTime); + + if (startWithFourier) + mode.SetDomain(EDomain::kFourier); + if (startWithAvg) + mode.SetAverage(true); + + return mode; +} + +//-------------------------------------------------------------------------- +// Prepare (public) +//-------------------------------------------------------------------------- +/** + *

Makes sure that everything the given mode needs has been calculated. + * + *

This is the single place which knows the calculation order and the + * dependencies between the representations. Note in particular that the + * difference of the Fourier spectra (f-d) and the Fourier transform of the + * difference (d-f) share the same histograms, distinguished by + * PMusrViewDataSet::diffFourierTag. Switching between the two therefore has to + * discard the previous result - something the former implementation did not do + * reliably, so that a d-f spectrum could survive into an f-d view. + * + * \param mode the view mode which shall be shown + * + * \return true if the mode can be shown + */ +Bool_t PMusrViewModel::Prepare(const PMusrViewMode &mode) +{ + if (!fValid) + return false; + + if (IsNonMusr()) { + if (mode.IsDifference() && !fDifferencePrepared) { + HandleDifference(); + fDifferencePrepared = true; + } + return true; + } + + if (fData.empty()) + return false; + + const EDiffStage stage = mode.GetDiffStage(); + + // the time domain difference is needed for the time difference view and as + // the input of the d-f Fourier transform + if (mode.IsDifference() && (stage == EDiffStage::kBeforeFourier)) { + if (!fDifferencePrepared || (fData[0].diff == nullptr)) { + HandleDifference(); + fDifferencePrepared = true; + } + } + + if (mode.IsFourier()) { + if (mode.IsDifference()) { + // d-f and f-d share diffFourier*, hence discard a stale one + const UInt_t wantedTag = (stage == EDiffStage::kBeforeFourier) ? 1 : 2; + if ((fData[0].diffFourierRe != nullptr) && (fData[0].diffFourierTag != wantedTag)) { + CleanupFourierDifference(); + fDiffFourierPrepared = false; + fFourierDiffPrepared = false; + fAveragePrepared = false; + } + + if (stage == EDiffStage::kBeforeFourier) { + if (!fDiffFourierPrepared || (fData[0].diffFourierRe == nullptr)) { + HandleDifferenceFourier(); + fDiffFourierPrepared = true; + } + } else { + if (!fFourierPrepared || (fData[0].dataFourierRe == nullptr)) { + HandleFourier(mode); + fFourierPrepared = true; + } + if (!fFourierDiffPrepared || (fData[0].diffFourierRe == nullptr)) { + HandleFourierDifference(); + fFourierDiffPrepared = true; + } + } + } else { + if (!fFourierPrepared || (fData[0].dataFourierRe == nullptr)) { + HandleFourier(mode); + fFourierPrepared = true; + } + } + + // the phase optimized real FT may be asked for later than the initial FT + if ((mode.GetFourierComponent() == EFourierComponent::kPhaseOptReal) && + (fData[0].dataFourierRe != nullptr) && + (fData[0].dataFourierPhaseOptReal == nullptr)) { + CalcPhaseOptReFT(); + fAveragePrepared = false; + } + } + + // the average is taken over whatever exists, hence it has to come last and + // has to be redone whenever the underlying representation changed + if (mode.IsAverage()) { + if (!fAveragePrepared || (fAverageKey != mode.GetKey()) || (fDataAvg.data == nullptr)) { + HandleAverage(); + fAveragePrepared = true; + fAverageKey = mode.GetKey(); + } + } + + return true; +} + +//-------------------------------------------------------------------------- +// SelectHisto (private) +//-------------------------------------------------------------------------- +/** + *

Picks the one histogram of a data set which the given mode shows. + * + *

This little function is what replaces the deeply nested mode cascades + * which used to be spelled out separately in the plotting routines and once + * more in the ascii export. + * + * \param dataSet the data set to pick from + * \param mode the view mode + * \param theory true to get the theory curve, false for the data curve + * \param imagPart true to get the imaginary part of the real+imaginary view + * + * \return the histogram, or nullptr if the mode shows no such curve + */ +TH1F *PMusrViewModel::SelectHisto(const PMusrViewDataSet &ds, const PMusrViewMode &mode, + Bool_t theory, Bool_t imagPart) const +{ + // in a difference view there is only the residual, no separate theory curve + if (mode.IsDifference()) { + if (theory) + return nullptr; + if (mode.IsTimeDomain()) + return imagPart ? nullptr : ds.diff; + switch (mode.GetFourierComponent()) { + case EFourierComponent::kReal: return imagPart ? nullptr : ds.diffFourierRe; + case EFourierComponent::kImag: return imagPart ? nullptr : ds.diffFourierIm; + case EFourierComponent::kRealAndImag: return imagPart ? ds.diffFourierIm : ds.diffFourierRe; + case EFourierComponent::kPower: return imagPart ? nullptr : ds.diffFourierPwr; + case EFourierComponent::kPhase: return imagPart ? nullptr : ds.diffFourierPhase; + case EFourierComponent::kPhaseOptReal: return imagPart ? nullptr : ds.diffFourierPhaseOptReal; + } + return nullptr; + } + + if (mode.IsTimeDomain()) { + if (imagPart) + return nullptr; + return theory ? ds.theory : ds.data; + } + + switch (mode.GetFourierComponent()) { + case EFourierComponent::kReal: + return imagPart ? nullptr : (theory ? ds.theoryFourierRe : ds.dataFourierRe); + case EFourierComponent::kImag: + return imagPart ? nullptr : (theory ? ds.theoryFourierIm : ds.dataFourierIm); + case EFourierComponent::kRealAndImag: + if (imagPart) + return theory ? ds.theoryFourierIm : ds.dataFourierIm; + return theory ? ds.theoryFourierRe : ds.dataFourierRe; + case EFourierComponent::kPower: + return imagPart ? nullptr : (theory ? ds.theoryFourierPwr : ds.dataFourierPwr); + case EFourierComponent::kPhase: + return imagPart ? nullptr : (theory ? ds.theoryFourierPhase : ds.dataFourierPhase); + case EFourierComponent::kPhaseOptReal: + return imagPart ? nullptr : (theory ? ds.theoryFourierPhaseOptReal : ds.dataFourierPhaseOptReal); + } + + return nullptr; +} + +//-------------------------------------------------------------------------- +// GetCurves (public) +//-------------------------------------------------------------------------- +/** + *

Collects all curves which the given mode makes visible. Call Prepare() + * first, otherwise the requested representation may not exist yet. + * + * \param mode the view mode + * + * \return the visible curves + */ +PMusrViewCurves PMusrViewModel::GetCurves(const PMusrViewMode &mode) const +{ + PMusrViewCurves curves; + + if (!fValid || IsNonMusr()) + return curves; + + // the y-errors are only meaningful for the time domain data + curves.fHasError = mode.IsTimeDomain(); + + if (mode.IsAverage()) { + TH1F *d = SelectHisto(fDataAvg, mode, kFALSE, kFALSE); + if (d == nullptr) + return curves; + curves.fData.push_back(d); + curves.fTheory.push_back(SelectHisto(fDataAvg, mode, kTRUE, kFALSE)); + curves.fDataSecond.push_back(SelectHisto(fDataAvg, mode, kFALSE, kTRUE)); + curves.fTheorySecond.push_back(SelectHisto(fDataAvg, mode, kTRUE, kTRUE)); + curves.fLabel.push_back("average"); + curves.fValid = true; + return curves; + } + + for (size_t i=0; iThe non-muSR error graphs which the given mode shows. Non-muSR data cannot + * be Fourier transformed, hence only the data and the difference view exist. + * + * \param mode the view mode + * \param theory true to get the theory graphs + * + * \return the visible error graphs + */ +std::vector PMusrViewModel::GetNonMusrCurves(const PMusrViewMode &mode, Bool_t theory) const +{ + std::vector result; + + if (!IsNonMusr()) + return result; + + for (const auto &ds : fNonMusrData) { + if (mode.IsDifference()) { + if (!theory && ds.diff) + result.push_back(ds.diff); + } else { + TGraphErrors *g = theory ? ds.theory : ds.data; + if (g) + result.push_back(g); + } + } + + return result; +} + +//-------------------------------------------------------------------------- +// GetFourierPhaseText (public) +//-------------------------------------------------------------------------- +/** + *

Describes the currently applied Fourier phase. + * + * \param mode the view mode + * + * \return the text, empty if no phase applies in this mode + */ +std::string PMusrViewModel::GetFourierPhaseText(const PMusrViewMode &mode) const +{ + // only the real/imaginary representations show the phase they were rotated by + const EFourierComponent comp = mode.GetFourierComponent(); + const Bool_t complexView = (comp == EFourierComponent::kReal) || + (comp == EFourierComponent::kImag) || + (comp == EFourierComponent::kRealAndImag); + if (!mode.IsFourier() || !complexView || fCurrentFourierPhase.empty()) + return std::string(); + + char buf[128]; + if (fCurrentFourierPhase.size() == 1) { + snprintf(buf, sizeof(buf), "phase = %.1f", fCurrentFourierPhase[0]); + return std::string(buf); + } + + std::string txt("phase = "); + for (size_t i=0; iThe x-range which the given mode shows by default, i.e. the msr-file range + * in the time domain and the Fourier plot range in the frequency domain. + * + * \param mode the view mode + * \param xmin where the minimum is returned + * \param xmax where the maximum is returned + */ +void PMusrViewModel::GetDefaultXRange(const PMusrViewMode &mode, Double_t &xmin, Double_t &xmax) const +{ + if (mode.IsTimeDomain()) { + xmin = fXmin; + xmax = fXmax; + return; + } + + if ((fFourier.fPlotRange[0] != -1.0) && (fFourier.fPlotRange[1] != -1.0)) { + xmin = fFourier.fPlotRange[0]; + xmax = fFourier.fPlotRange[1]; + return; + } + + // fall back onto the full range of the transform + xmin = 0.0; + xmax = 0.0; + const PMusrViewCurves curves = GetCurves(mode); + if (!curves.IsEmpty() && curves.fData[0]) { + TH1F *h = curves.fData[0]; + xmin = h->GetBinLowEdge(1); + xmax = h->GetBinLowEdge(h->GetNbinsX()) + h->GetBinWidth(h->GetNbinsX()); + } +} + +//-------------------------------------------------------------------------- +// GetDefaultYRange (public) +//-------------------------------------------------------------------------- +/** + *

The y-range which the given mode shows by default for the given x-range. + * An explicit y-range of the msr-file PLOT block wins in the time domain. + * + * \param mode the view mode + * \param xmin lower bound of the x-range + * \param xmax upper bound of the x-range + * \param ymin where the minimum is returned + * \param ymax where the maximum is returned + */ +void PMusrViewModel::GetDefaultYRange(const PMusrViewMode &mode, Double_t xmin, Double_t xmax, + Double_t &ymin, Double_t &ymax) const +{ + ymin = 0.0; + ymax = 0.0; + + if (mode.IsTimeDomain() && !mode.IsDifference() && fYRangePresent) { + ymin = fYmin; + ymax = fYmax; + return; + } + + PMusrViewModel *self = const_cast(this); + + // An averaged spectrum must not rescale the axis: the average of several + // data sets is meant to be compared with them, hence the range is the one + // of the individual spectra of the very same representation. + PMusrViewMode rangeMode(mode); + rangeMode.SetAverage(false); + + std::vector histos; + if (fValid && !IsNonMusr()) { + const PMusrViewCurves curves = GetCurves(rangeMode); + for (auto *h : curves.fData) + if (h) histos.push_back(h); + for (auto *h : curves.fDataSecond) + if (h) histos.push_back(h); + // the theory has to fit into the frame as well + for (auto *h : curves.fTheory) + if (h) histos.push_back(h); + for (auto *h : curves.fTheorySecond) + if (h) histos.push_back(h); + } + + if (histos.empty()) { + // BuildData() calls this before any curve exists + for (const auto &ds : fData) + if (ds.data) histos.push_back(ds.data); + } + + if (histos.empty()) + return; + + // The two representations use different conventions, which are kept as they + // were: in the time domain the extrema are taken over the visible x-range and + // widened by an absolute margin, in the frequency domain they are taken over + // the complete spectrum and widened by a relative one. + if (rangeMode.IsFourier()) { + ymin = self->GetMinimum(histos[0]); + ymax = self->GetMaximum(histos[0]); + for (size_t i=1; iGetMinimum(histos[i]); + const Double_t hi = self->GetMaximum(histos[i]); + if (lo < ymin) ymin = lo; + if (hi > ymax) ymax = hi; + } + // the power spectrum is positive, hence its lower bound is pulled down + // rather than pushed further away from zero + ymin *= (rangeMode.GetFourierComponent() == EFourierComponent::kPower) ? 0.95 : 1.05; + ymax *= 1.05; + } else { + ymin = self->GetMinimum(histos[0], xmin, xmax); + ymax = self->GetMaximum(histos[0], xmin, xmax); + for (size_t i=1; iGetMinimum(histos[i], xmin, xmax); + const Double_t hi = self->GetMaximum(histos[i], xmin, xmax); + if (lo < ymin) ymin = lo; + if (hi > ymax) ymax = hi; + } + + const Double_t dd = 0.05*fabs(ymax-ymin); + ymin -= dd; + ymax += dd; + } + + if (fMsrHandler && fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fLogY && (ymax > 0.0)) + ymin = 1.0e-4 * ymax; +} + +//-------------------------------------------------------------------------- +// GetAxisTitles (public) +//-------------------------------------------------------------------------- +/** + *

The axis titles belonging to the given mode. + * + * \param mode the view mode + * \param xTitle where the x-axis title is returned + * \param yTitle where the y-axis title is returned + */ +void PMusrViewModel::GetAxisTitles(const PMusrViewMode &mode, std::string &xTitle, std::string &yTitle) const +{ + const Bool_t avg = mode.IsAverage(); + + if (IsNonMusr()) { + // non-muSR data carry their own axis labels, which the data handler took + // from the data file; there is no time or frequency axis here + xTitle.clear(); + yTitle.clear(); + + if (fMsrHandler && fRunList) { + PMsrPlotStructure plotInfo = fMsrHandler->GetMsrPlotList()->at(fPlotNumber); + if (!plotInfo.fRuns.empty()) { + const UInt_t runNo = static_cast(plotInfo.fRuns[0])-1; + const TString runName = *fMsrHandler->GetMsrRunList()->at(runNo).GetRunName(); + xTitle = fRunList->GetXAxisTitle(runName, runNo); + yTitle = fRunList->GetYAxisTitle(runName, runNo); + } + } + + if (mode.IsDifference()) + yTitle = avg ? "" : "data-theory"; + + return; + } + + if (mode.IsTimeDomain()) { + if (fPlotType == MSR_PLOT_BNMR) + xTitle = "time (s)"; // BNMR/BNQR runs use seconds + else + xTitle = "time (#mus)"; + + if (mode.IsDifference()) { + yTitle = avg ? "" : "data-theory"; + return; + } + + switch (fPlotType) { + case MSR_PLOT_SINGLE_HISTO: + if (fMsrHandler && fMsrHandler->GetMsrRunList()->at(0).IsLifetimeCorrected()) + yTitle = avg ? "" : "Asymmetry"; + else if (fScaleN0AndBkg) + yTitle = avg ? " per nsec" : "N(t) per nsec"; + else + yTitle = avg ? " per bin" : "N(t) per bin"; + break; + case MSR_PLOT_SINGLE_HISTO_RRF: + case MSR_PLOT_ASYM_RRF: + yTitle = avg ? "" : "RRF Asymmetry"; + break; + case MSR_PLOT_MU_MINUS: + yTitle = avg ? " per bin" : "N(t) per bin"; + break; + case MSR_PLOT_BNMR: + case MSR_PLOT_ASYM: + yTitle = avg ? "" : "Asymmetry"; + break; + default: + yTitle = "??"; + break; + } + return; + } + + // frequency domain + switch (fFourier.fUnits) { + case FOURIER_UNIT_GAUSS: xTitle = "Field (G)"; break; + case FOURIER_UNIT_TESLA: xTitle = "Field (T)"; break; + case FOURIER_UNIT_FREQ: xTitle = "Frequency (MHz)"; break; + case FOURIER_UNIT_CYCLES: xTitle = "Frequency (Mc/s)";break; + default: xTitle = "??"; break; + } + + switch (mode.GetFourierComponent()) { + case EFourierComponent::kReal: yTitle = "Real Fourier"; break; + case EFourierComponent::kImag: yTitle = "Imaginary Fourier"; break; + case EFourierComponent::kRealAndImag: yTitle = "Real+Imag Fourier"; break; + case EFourierComponent::kPower: yTitle = "Ampl. Fourier"; break; + case EFourierComponent::kPhase: yTitle = "Phase Fourier"; break; + case EFourierComponent::kPhaseOptReal: yTitle = "Phase Opt. Real Fourier"; break; + } + + // say which of the two difference flavours is shown + switch (mode.GetDiffStage()) { + case EDiffStage::kNone: + break; + case EDiffStage::kBeforeFourier: + yTitle += " (d-f: data-theory)"; + break; + case EDiffStage::kAfterFourier: + yTitle += " (f-d: [F data]-[F theory])"; + break; + } + + if (avg) + yTitle = "<" + yTitle + ">"; +} + +//-------------------------------------------------------------------------- +// IsLogX (public) +//-------------------------------------------------------------------------- +/** + * \return true if the msr-file PLOT block asks for a logarithmic x-axis + */ +Bool_t PMusrViewModel::IsLogX() const +{ + if (fMsrHandler == nullptr) + return kFALSE; + return fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fLogX; +} + +//-------------------------------------------------------------------------- +// IsLogY (public) +//-------------------------------------------------------------------------- +/** + * \return true if the msr-file PLOT block asks for a logarithmic y-axis + */ +Bool_t PMusrViewModel::IsLogY() const +{ + if (fMsrHandler == nullptr) + return kFALSE; + return fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fLogY; +} + +//-------------------------------------------------------------------------- +// BuildTextBlocks (private) +//-------------------------------------------------------------------------- +/** + *

Generates the title, the fit parameter block, the theory/functions block + * and the run info block as plain text. + * + *

The former implementation wrote these directly into TPaveText/TLegend + * objects, which scaled the font down until the content became unreadable. + * Producing plain text here lets the front end decide how to present it - the + * web front end puts each block into its own scrollable panel. + */ +void PMusrViewModel::BuildTextBlocks() +{ + fParameterLines.clear(); + fTheoryLines.clear(); + fInfoLines.clear(); + + if (fMsrHandler == nullptr) + return; + + fTitle = fMsrHandler->GetMsrTitle()->Data(); + + BuildRRFText(); + + //-- fit parameters ------------------------------------------------------ + PMsrParamList param = *fMsrHandler->GetMsrParamList(); + + Int_t maxLength = 0; + for (UInt_t i=0; i maxLength) + maxLength = param[i].fName.Length(); + } + maxLength += 2; + + Char_t cnum[128]; + Char_t accStr[32]; + for (UInt_t i=0; i(strlen(cnum)); j++) + str += " "; + str += " "; + // parameter error + if (param[i].fPosErrorPresent) { // minos was used + Double_t err = (param[i].fPosError - param[i].fStep) / 2.0; + if ((fabs(fabs(param[i].fStep) - param[i].fPosError) < 0.1*fabs(param[i].fStep)) && + (fabs(fabs(param[i].fStep) - param[i].fPosError) < 0.1*param[i].fPosError)) { + if (round(err)-err==0) + snprintf(cnum, sizeof(cnum), "%.1lf", err); + else + snprintf(cnum, sizeof(cnum), accStr, err); + } else { + snprintf(accStr, sizeof(accStr), "%%.%dlf!!", accuracy); + if (round(err)-err==0) + snprintf(cnum, sizeof(cnum), "%.1lf!!", err); + else + snprintf(cnum, sizeof(cnum), accStr, err); + } + str += cnum; + } else { // minos was not used + if (round(param[i].fStep)-param[i].fStep==0) + snprintf(cnum, sizeof(cnum), "%.1lf", param[i].fStep); + else + snprintf(cnum, sizeof(cnum), accStr, param[i].fStep); + str += cnum; + } + fParameterLines.push_back(str.Data()); + } + + //-- theory -------------------------------------------------------------- + PMsrLines theory = *fMsrHandler->GetMsrTheory(); + for (UInt_t i=1; i 0) { // strip the comment + str.Resize(idx-1); + str.Resize(str.Strip().Length()); + } + fTheoryLines.push_back(str.Data()); + } + + //-- functions ----------------------------------------------------------- + PMsrLines functions = *fMsrHandler->GetMsrFunctions(); + if (functions.size() > 1) + fTheoryLines.push_back(""); + for (UInt_t i=1; iGetMsrStatistic(); + Char_t num[64]; + + TString tstr = "musrfit: "; + tstr += statistic.fDate; + tstr += TString(", "); + tstr += statistic.fChisq ? TString("chisq = ") : TString("maxLH = "); + snprintf(num, sizeof(num), "%.10g", statistic.fMin); + tstr += num; + tstr += TString(" , NDF = "); + tstr += statistic.fNdf; + tstr += statistic.fChisq ? TString(" , chisq/NDF = ") : TString(" , maxLH/NDF = "); + if (statistic.fNdf != 0) { + snprintf(num, sizeof(num), "%.5g", statistic.fMin/statistic.fNdf); + tstr += num; + } else { + tstr += TString("undefined"); + } + fStatisticLine = tstr.Data(); + + //-- run info, one line per data set ------------------------------------ + if (IsNonMusr()) { + // for non-muSR data the runs are identified by what is plotted against + // what, exactly as the legend of the former implementation did + PMsrPlotStructure plotInfoNonMusr = fMsrHandler->GetMsrPlotList()->at(fPlotNumber); + for (UInt_t i=0; i(plotInfoNonMusr.fRuns[i])-1; + const TString runName = *fMsrHandler->GetMsrRunList()->at(runNo).GetRunName(); + TString label = fRunList->GetYAxisTitle(runName, runNo); + label += " vs. "; + label += fRunList->GetXAxisTitle(runName, runNo); + fInfoLines.push_back(label.Data()); + } + return; + } + + double dval; + const PDoublePairVector *ddvec; + Char_t sval[128]; + PMsrPlotStructure plotInfo = fMsrHandler->GetMsrPlotList()->at(fPlotNumber); + PMsrRunList runs = *fMsrHandler->GetMsrRunList(); + + for (UInt_t i=0; i(plotInfo.fRuns[i])-1; + if (runs[runNo].GetRunNameSize() > 1) + tstr = "++" + *runs[runNo].GetRunName() + TString(","); + else + tstr = *runs[runNo].GetRunName() + TString(","); + + 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) || + (runs[runNo].GetFitType() == MSR_FITTYPE_ASYM_RRF) || + (runs[runNo].GetFitType() == MSR_FITTYPE_BNMR)) { + tstr += TString("h:"); + TString grouping; + fMsrHandler->GetGroupingString(runNo, "forward", grouping); + tstr += grouping; + tstr += TString("/"); + grouping = ""; + fMsrHandler->GetGroupingString(runNo, "backward", grouping); + tstr += grouping; + tstr += TString(","); + } + + ddvec = fRunList->GetTemp(*runs[runNo].GetRunName()); + if (ddvec->empty()) { + tstr += TString("T=??,"); + } else if (ddvec->size() == 1) { + tstr += TString("T="); + snprintf(sval, sizeof(sval), "%0.2lf", ddvec->at(0).first); + tstr += TString(sval) + TString("K,"); + } else { + for (UInt_t j=0; jsize(); ++j) { + snprintf(sval, sizeof(sval), "T%u=", j); + tstr += TString(sval); + snprintf(sval, sizeof(sval), "%0.2lf", ddvec->at(j).first); + tstr += TString(sval) + TString("K,"); + } + } + + tstr += TString("B="); + dval = fRunList->GetField(*runs[runNo].GetRunName()); + if (dval == PMUSR_UNDEFINED) { + tstr += TString("??,"); + } else if (dval < 1.0e4) { // Gauss makes sense as a unit + snprintf(sval, sizeof(sval), "%0.2lf", dval); + tstr += TString(sval) + TString("G,"); + } else { // Tesla makes sense as a unit + snprintf(sval, sizeof(sval), "%0.2lf", dval/1.0e4); + tstr += TString(sval) + TString("T,"); + } + + tstr += TString("E="); + dval = fRunList->GetEnergy(*runs[runNo].GetRunName()); + if (dval == PMUSR_UNDEFINED) { + tstr += TString("??,"); + } else if (dval < 1.0e3) { // keV makes sense as a unit + snprintf(sval, sizeof(sval), "%0.2lf", dval); + tstr += TString(sval) + TString("keV,"); + } else { // MeV makes sense as a unit + snprintf(sval, sizeof(sval), "%0.2lf", dval/1.0e3); + tstr += TString(sval) + TString("MeV,"); + } + + tstr += fRunList->GetSetup(*runs[runNo].GetRunName()); + + fInfoLines.push_back(tstr.Data()); + } +} + +//-------------------------------------------------------------------------- +// ExportData (public) +//-------------------------------------------------------------------------- +/** + *

Writes an ascii dump of exactly what the given mode shows, restricted to + * the x-range [xmin, xmax]. + * + *

The former implementation repeated the complete view mode cascade here, + * some 700 lines duplicating the plotting logic. Going through GetCurves() + * means the export can never again disagree with what is on screen. + * + * \param fileName name of the ascii file to be written + * \param mode the view mode + * \param xmin lower bound of the range to be exported + * \param xmax upper bound of the range to be exported + * + * \return true on success + */ +Bool_t PMusrViewModel::ExportData(const char *fileName, const PMusrViewMode &mode, + Double_t xmin, Double_t xmax) +{ + if (fileName == nullptr) { + fLastError = "no path file name provided."; + std::cerr << std::endl << ">> PMusrViewModel::ExportData(): **ERROR** " << fLastError << std::endl; + return false; + } + + PMusrViewAsciiDumpVector dumpVector; + std::vector header; + + //-- column naming, derived from the view mode; the names are the ones the + //-- former hand written export produced, so that existing scripts keep + //-- working -------------------------------------------------------------- + std::string xName, xTheoName, valueName, errName, theoName, secondName, secondErrName; + + if (mode.IsFourier()) { + std::string comp; + switch (mode.GetFourierComponent()) { + case EFourierComponent::kReal: comp = "Re"; break; + case EFourierComponent::kImag: comp = "Im"; break; + case EFourierComponent::kRealAndImag: comp = "Re"; break; + case EFourierComponent::kPower: comp = "Pwr"; break; + case EFourierComponent::kPhase: comp = "Phase"; break; + case EFourierComponent::kPhaseOptReal: comp = "PhaseOptRe"; break; + } + xName = mode.IsDifference() ? "x" : "freq"; + xTheoName = "freqTheo"; + valueName = (mode.IsDifference() ? "F_diff" : "F_") + comp; + errName = "err" + valueName; + theoName = "F_theo" + comp; + secondName = mode.IsDifference() ? "F_diffIm" : "F_Im"; + secondErrName = "err" + secondName; + } else { + xName = mode.IsDifference() ? "x" : "xData"; + xTheoName = "xTheory"; + valueName = mode.IsDifference() ? "diff" : "data"; + errName = mode.IsDifference() ? "errDiff" : "errData"; + theoName = "theory"; + secondName = valueName; + secondErrName = errName; + } + + //-- collect the columns ------------------------------------------------- + if (IsNonMusr()) { + std::vector data = GetNonMusrCurves(mode, kFALSE); + std::vector theo = GetNonMusrCurves(mode, kTRUE); + + for (size_t i=0; iGetN(); j++) { + data[i]->GetPoint(j, x, y); + if ((x < xmin) || (x > xmax)) + continue; + dump.dataX.push_back(x); + dump.data.push_back(y); + dump.dataErr.push_back(data[i]->GetErrorY(j)); + } + dumpVector.push_back(dump); + header.push_back("xData" + std::to_string(i)); + header.push_back("data" + std::to_string(i)); + header.push_back("errData" + std::to_string(i)); + + if (i < theo.size() && theo[i]) { + PMusrViewAsciiDump th; + for (Int_t j=0; jGetN(); j++) { + theo[i]->GetPoint(j, x, y); + if ((x < xmin) || (x > xmax)) + continue; + th.dataX.push_back(x); + th.data.push_back(y); + } + dumpVector.push_back(th); + header.push_back("xTheory" + std::to_string(i)); + header.push_back("theory" + std::to_string(i)); + } + } + } else { + const PMusrViewCurves curves = GetCurves(mode); + if (!curves.fValid) { + fLastError = "nothing to export for the current view."; + std::cerr << std::endl << ">> PMusrViewModel::ExportData(): **ERROR** " << fLastError << std::endl; + return false; + } + + // a data column is always followed by its theory column, as before + for (size_t i=0; i> PMusrViewModel::ExportData(): **ERROR** " << fLastError << std::endl; + return false; + } + + std::ofstream fout(fileName, std::ofstream::out); + if (!fout.is_open()) { + fLastError = std::string("couldn't open '") + fileName + "' for writing."; + std::cerr << std::endl << ">> PMusrViewModel::ExportData(): **ERROR** " << fLastError << std::endl; + return false; + } + + //-- header -------------------------------------------------------------- + fout << "% " << fTitle << std::endl; + fout << "% view: " << mode.GetLabel() << std::endl; + if (mode.IsAverage()) + fout << "% from averaged view" << std::endl; + fout << "% "; + for (size_t i=0; i maxLength) + maxLength = d.dataX.size(); + + for (size_t row=0; row> PMusrViewModel::ExportData(): data exported to '" << fileName << "'." << std::endl; + + return true; +} + + +//-------------------------------------------------------------------------- +// SetMsrHandler (public) +//-------------------------------------------------------------------------- +/** + *

Keeps the msr-handler and takes over the settings of the FOURIER block. + * + * \param msrHandler pointer to the msr-file handler + */ +void PMusrViewModel::SetMsrHandler(PMsrHandler *msrHandler) +{ + fMsrHandler = msrHandler; + + fScaleN0AndBkg = IsScaleN0AndBkg(); + + if (fMsrHandler->GetMsrFourierList()->fFourierBlockPresent) { + fFourier.fFourierBlockPresent = true; + if (fMsrHandler->GetMsrFourierList()->fUnits != FOURIER_UNIT_NOT_GIVEN) + fFourier.fUnits = fMsrHandler->GetMsrFourierList()->fUnits; + if (fMsrHandler->GetMsrFourierList()->fFourierPower != -1) + fFourier.fFourierPower = fMsrHandler->GetMsrFourierList()->fFourierPower; + fFourier.fDCCorrected = fMsrHandler->GetMsrFourierList()->fDCCorrected; + if (fMsrHandler->GetMsrFourierList()->fApodization != FOURIER_APOD_NOT_GIVEN) + fFourier.fApodization = fMsrHandler->GetMsrFourierList()->fApodization; + if (fMsrHandler->GetMsrFourierList()->fPlotTag != FOURIER_PLOT_NOT_GIVEN) + fFourier.fPlotTag = fMsrHandler->GetMsrFourierList()->fPlotTag; + fFourier.fPhase = fMsrHandler->GetMsrFourierList()->fPhase; + + if ((fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0] != -1.0) && + (fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1] != -1.0)) { + fFourier.fRangeForPhaseCorrection[0] = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; + fFourier.fRangeForPhaseCorrection[1] = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; + } + if ((fMsrHandler->GetMsrFourierList()->fPlotRange[0] != -1.0) && + (fMsrHandler->GetMsrFourierList()->fPlotRange[1] != -1.0)) { + fFourier.fPlotRange[0] = fMsrHandler->GetMsrFourierList()->fPlotRange[0]; + fFourier.fPlotRange[1] = fMsrHandler->GetMsrFourierList()->fPlotRange[1]; + } + } +} + +//-------------------------------------------------------------------------- +// BuildRRFText (private) +//-------------------------------------------------------------------------- +/** + *

Assembles the RRF information string, if the msr-file uses a rotating + * reference frame. + */ +void PMusrViewModel::BuildRRFText() +{ + fRRFText.clear(); + + if (fMsrHandler == nullptr) + return; + + const Bool_t rrfPlot = (fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) && + (fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0); + const Bool_t rrfFit = (fMsrHandler->GetMsrGlobal()->GetRRFPacking() > 0) && + fMsrHandler->GetMsrGlobal()->GetRRFUnit().CompareTo("??"); + + if (!rrfPlot && !rrfFit) + return; + + Int_t rrfUnitTag = -1; + Double_t rrfFreq = 0.0; + Int_t rrfPacking = 0; + + if (rrfPlot) { // RRF single histo PLOT + rrfUnitTag = fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit; + rrfFreq = fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + rrfPacking = fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking; + } else { // RRF single histo FIT + rrfUnitTag = fMsrHandler->GetMsrGlobal()->GetRRFUnitTag(); + rrfFreq = fMsrHandler->GetMsrGlobal()->GetRRFFreq(fMsrHandler->GetMsrGlobal()->GetRRFUnit().Data()); + rrfPacking = fMsrHandler->GetMsrGlobal()->GetRRFPacking(); + } + + Char_t freqStr[64]; + snprintf(freqStr, sizeof(freqStr), "%.5g", rrfFreq); + + TString txt("RRF: "); + switch (rrfUnitTag) { + case RRF_UNIT_kHz: txt += TString("#nu_{RRF} = ") + freqStr + TString(" (kHz)"); break; + case RRF_UNIT_MHz: txt += TString("#nu_{RRF} = ") + freqStr + TString(" (MHz)"); break; + case RRF_UNIT_Mcs: txt += TString("#omega_{RRF} = ") + freqStr + TString(" (Mc/s)"); break; + case RRF_UNIT_G: txt += TString("B_{RRF} = ") + freqStr + TString(" (G)"); break; + case RRF_UNIT_T: txt += TString("B_{RRF} = ") + freqStr + TString(" (T)"); break; + default: break; + } + txt += TString(", RRF packing = "); + txt += rrfPacking; + + fRRFText = txt.Data(); +} + +//-------------------------------------------------------------------------- +// IncrementFourierPhase (public) +//-------------------------------------------------------------------------- +/** + *

Rotates the complex Fourier spectra of data, theory and difference by the + * configured phase increment. + * + * \param mode the current view mode, used to decide whether a phase applies + */ +void PMusrViewModel::IncrementFourierPhase(const PMusrViewMode &mode) +{ + if (!mode.IsPhaseAdjustable()) + return; + + Double_t re, im; + const Double_t cp = TMath::Cos(fFourier.fPhaseIncrement/180.0*TMath::Pi()); + const Double_t sp = TMath::Sin(fFourier.fPhaseIncrement/180.0*TMath::Pi()); + + for (UInt_t i=0; iGetNbinsX(); j++) { + re = fData[i].dataFourierRe->GetBinContent(j) * cp + fData[i].dataFourierIm->GetBinContent(j) * sp; + im = fData[i].dataFourierIm->GetBinContent(j) * cp - fData[i].dataFourierRe->GetBinContent(j) * sp; + fData[i].dataFourierRe->SetBinContent(j, re); + fData[i].dataFourierIm->SetBinContent(j, im); + } + } + if ((fData[i].theoryFourierRe != nullptr) && (fData[i].theoryFourierIm != nullptr)) { + for (Int_t j=0; jGetNbinsX(); j++) { + re = fData[i].theoryFourierRe->GetBinContent(j) * cp + fData[i].theoryFourierIm->GetBinContent(j) * sp; + im = fData[i].theoryFourierIm->GetBinContent(j) * cp - fData[i].theoryFourierRe->GetBinContent(j) * sp; + fData[i].theoryFourierRe->SetBinContent(j, re); + fData[i].theoryFourierIm->SetBinContent(j, im); + } + } + if ((fData[i].diffFourierRe != nullptr) && (fData[i].diffFourierIm != nullptr)) { + for (Int_t j=0; jGetNbinsX(); j++) { + re = fData[i].diffFourierRe->GetBinContent(j) * cp + fData[i].diffFourierIm->GetBinContent(j) * sp; + im = fData[i].diffFourierIm->GetBinContent(j) * cp - fData[i].diffFourierRe->GetBinContent(j) * sp; + fData[i].diffFourierRe->SetBinContent(j, re); + fData[i].diffFourierIm->SetBinContent(j, im); + } + } + } + + // the average was taken over the unrotated spectra + fAveragePrepared = false; +} + +//-------------------------------------------------------------------------- +// DecrementFourierPhase (public) +//-------------------------------------------------------------------------- +/** + *

Rotates the complex Fourier spectra of data, theory and difference by the + * negative configured phase increment. + * + * \param mode the current view mode, used to decide whether a phase applies + */ +void PMusrViewModel::DecrementFourierPhase(const PMusrViewMode &mode) +{ + if (!mode.IsPhaseAdjustable()) + return; + + Double_t re, im; + const Double_t cp = TMath::Cos(fFourier.fPhaseIncrement/180.0*TMath::Pi()); + const Double_t sp = TMath::Sin(fFourier.fPhaseIncrement/180.0*TMath::Pi()); + + for (UInt_t i=0; iGetNbinsX(); j++) { + re = fData[i].dataFourierRe->GetBinContent(j) * cp - fData[i].dataFourierIm->GetBinContent(j) * sp; + im = fData[i].dataFourierIm->GetBinContent(j) * cp + fData[i].dataFourierRe->GetBinContent(j) * sp; + fData[i].dataFourierRe->SetBinContent(j, re); + fData[i].dataFourierIm->SetBinContent(j, im); + } + } + if ((fData[i].theoryFourierRe != nullptr) && (fData[i].theoryFourierIm != nullptr)) { + for (Int_t j=0; jGetNbinsX(); j++) { + re = fData[i].theoryFourierRe->GetBinContent(j) * cp - fData[i].theoryFourierIm->GetBinContent(j) * sp; + im = fData[i].theoryFourierIm->GetBinContent(j) * cp + fData[i].theoryFourierRe->GetBinContent(j) * sp; + fData[i].theoryFourierRe->SetBinContent(j, re); + fData[i].theoryFourierIm->SetBinContent(j, im); + } + } + if ((fData[i].diffFourierRe != nullptr) && (fData[i].diffFourierIm != nullptr)) { + for (Int_t j=0; jGetNbinsX(); j++) { + re = fData[i].diffFourierRe->GetBinContent(j) * cp - fData[i].diffFourierIm->GetBinContent(j) * sp; + im = fData[i].diffFourierIm->GetBinContent(j) * cp + fData[i].diffFourierRe->GetBinContent(j) * sp; + fData[i].diffFourierRe->SetBinContent(j, re); + fData[i].diffFourierIm->SetBinContent(j, im); + } + } + } + + fAveragePrepared = false; +} + +//-------------------------------------------------------------------------- +// InitFourier (private) +//-------------------------------------------------------------------------- +/** + *

Initializes the Fourier structure. + */ +void PMusrViewModel::InitFourier() +{ + fFourier.fFourierBlockPresent = false; // fourier block present + fFourier.fUnits = FOURIER_UNIT_GAUSS; // fourier untis + fFourier.fFourierPower = 0; // no zero padding + fFourier.fApodization = FOURIER_APOD_NONE; // no apodization + fFourier.fPlotTag = FOURIER_PLOT_REAL_AND_IMAG; // initial plot tag, plot real and imaginary part + fFourier.fPhaseParamNo.clear(); + fFourier.fPhase.clear(); + for (UInt_t i=0; i<2; i++) { + fFourier.fRangeForPhaseCorrection[i] = -1.0; // frequency range for phase correction, default: {-1, -1} = NOT GIVEN + fFourier.fPlotRange[i] = -1.0; // fourier plot range, default: {-1, -1} = NOT GIVEN + } + fFourier.fPhaseIncrement = 1.0; // fourier phase increment +} + +//-------------------------------------------------------------------------- +// InitAverage (private) +//-------------------------------------------------------------------------- +/** + *

Initializes the Average structure. + */ +void PMusrViewModel::InitAverage() +{ + fDataAvg.data = nullptr; + fDataAvg.dataFourierRe = nullptr; + fDataAvg.dataFourierIm = nullptr; + fDataAvg.dataFourierPwr = nullptr; + fDataAvg.dataFourierPhase = nullptr; + fDataAvg.dataFourierPhaseOptReal = nullptr; + fDataAvg.theory = nullptr; + fDataAvg.theoryFourierRe = nullptr; + fDataAvg.theoryFourierIm = nullptr; + fDataAvg.theoryFourierPwr = nullptr; + fDataAvg.theoryFourierPhase = nullptr; + fDataAvg.theoryFourierPhaseOptReal = nullptr; + fDataAvg.diff = nullptr; + fDataAvg.diffFourierRe = nullptr; + fDataAvg.diffFourierIm = nullptr; + fDataAvg.diffFourierPwr = nullptr; + fDataAvg.diffFourierPhase = nullptr; + fDataAvg.diffFourierPhaseOptReal = nullptr; + fDataAvg.dataRange = nullptr; + fDataAvg.diffFourierTag = 0; +} + +//-------------------------------------------------------------------------- +// InitDataSet (private) +//-------------------------------------------------------------------------- +/** + *

Initializes the data set (histogram representation). + * + * \param dataSet data set to be initialized + */ +void PMusrViewModel::InitDataSet(PMusrViewDataSet &dataSet) +{ + dataSet.data = nullptr; + dataSet.dataFourierRe = nullptr; + dataSet.dataFourierIm = nullptr; + dataSet.dataFourierPwr = nullptr; + dataSet.dataFourierPhase = nullptr; + dataSet.dataFourierPhaseOptReal = nullptr; + dataSet.theory = nullptr; + dataSet.theoryFourierRe = nullptr; + dataSet.theoryFourierIm = nullptr; + dataSet.theoryFourierPwr = nullptr; + dataSet.theoryFourierPhase = nullptr; + dataSet.theoryFourierPhaseOptReal = nullptr; + dataSet.diff = nullptr; + dataSet.diffFourierRe = nullptr; + dataSet.diffFourierIm = nullptr; + dataSet.diffFourierPwr = nullptr; + dataSet.diffFourierPhase = nullptr; + dataSet.diffFourierPhaseOptReal = nullptr; + dataSet.dataRange = nullptr; +} + +//-------------------------------------------------------------------------- +// InitDataSet (private) +//-------------------------------------------------------------------------- +/** + *

Initializes the data set (error graph representation). + * + * \param dataSet data set to be initialized + */ +void PMusrViewModel::InitDataSet(PMusrViewNonMusrDataSet &dataSet) +{ + dataSet.data = nullptr; + dataSet.dataFourierRe = nullptr; + dataSet.dataFourierIm = nullptr; + dataSet.dataFourierPwr = nullptr; + dataSet.dataFourierPhase = nullptr; + dataSet.theory = nullptr; + dataSet.theoryFourierRe = nullptr; + dataSet.theoryFourierIm = nullptr; + dataSet.theoryFourierPwr = nullptr; + dataSet.theoryFourierPhase = nullptr; + dataSet.diff = nullptr; + dataSet.diffFourierRe = nullptr; + dataSet.diffFourierIm = nullptr; + dataSet.diffFourierPwr = nullptr; + dataSet.diffFourierPhase = nullptr; + dataSet.dataRange = nullptr; +} + +//-------------------------------------------------------------------------- +// CleanupDataSet (private) +//-------------------------------------------------------------------------- +/** + *

Cleans up a data set (histogram representation). + * + * \param dataSet data set to be cleaned up. + */ +void PMusrViewModel::CleanupDataSet(PMusrViewDataSet &dataSet) +{ + if (dataSet.data) { + delete dataSet.data; + dataSet.data = nullptr; + } + if (dataSet.dataFourierRe) { + delete dataSet.dataFourierRe; + dataSet.dataFourierRe = nullptr; + } + if (dataSet.dataFourierIm) { + delete dataSet.dataFourierIm; + dataSet.dataFourierIm = nullptr; + } + if (dataSet.dataFourierPwr) { + delete dataSet.dataFourierPwr; + dataSet.dataFourierPwr = nullptr; + } + if (dataSet.dataFourierPhase) { + delete dataSet.dataFourierPhase; + dataSet.dataFourierPhase = nullptr; + } + if (dataSet.dataFourierPhaseOptReal) { + delete dataSet.dataFourierPhaseOptReal; + dataSet.dataFourierPhaseOptReal = nullptr; + } + if (dataSet.theory) { + delete dataSet.theory; + dataSet.theory = nullptr; + } + if (dataSet.theoryFourierRe) { + delete dataSet.theoryFourierRe; + dataSet.theoryFourierRe = nullptr; + } + if (dataSet.theoryFourierIm) { + delete dataSet.theoryFourierIm; + dataSet.theoryFourierIm = nullptr; + } + if (dataSet.theoryFourierPwr) { + delete dataSet.theoryFourierPwr; + dataSet.theoryFourierPwr = nullptr; + } + if (dataSet.theoryFourierPhase) { + delete dataSet.theoryFourierPhase; + dataSet.theoryFourierPhase = nullptr; + } + if (dataSet.theoryFourierPhaseOptReal) { + delete dataSet.theoryFourierPhaseOptReal; + dataSet.theoryFourierPhaseOptReal = nullptr; + } + if (dataSet.diff) { + delete dataSet.diff; + dataSet.diff = nullptr; + } + if (dataSet.diffFourierRe) { + delete dataSet.diffFourierRe; + dataSet.diffFourierRe = nullptr; + } + if (dataSet.diffFourierIm) { + delete dataSet.diffFourierIm; + dataSet.diffFourierIm = nullptr; + } + if (dataSet.diffFourierPwr) { + delete dataSet.diffFourierPwr; + dataSet.diffFourierPwr = nullptr; + } + if (dataSet.diffFourierPhase) { + delete dataSet.diffFourierPhase; + dataSet.diffFourierPhase = nullptr; + } + if (dataSet.diffFourierPhaseOptReal) { + delete dataSet.diffFourierPhaseOptReal; + dataSet.diffFourierPhaseOptReal = nullptr; + } + if (dataSet.dataRange) { + delete dataSet.dataRange; + dataSet.dataRange = nullptr; + } +} + +//-------------------------------------------------------------------------- +// CleanupDataSet (private) +//-------------------------------------------------------------------------- +/** + *

Cleans up a data set (error graph representation). + * + * \param dataSet data set to be cleaned up. + */ +void PMusrViewModel::CleanupDataSet(PMusrViewNonMusrDataSet &dataSet) +{ + if (dataSet.data) { + delete dataSet.data; + dataSet.data = nullptr; + } + if (dataSet.dataFourierRe) { + delete dataSet.dataFourierRe; + dataSet.dataFourierRe = nullptr; + } + if (dataSet.dataFourierIm) { + delete dataSet.dataFourierIm; + dataSet.dataFourierIm = nullptr; + } + if (dataSet.dataFourierPwr) { + delete dataSet.dataFourierPwr; + dataSet.dataFourierPwr = nullptr; + } + if (dataSet.dataFourierPhase) { + delete dataSet.dataFourierPhase; + dataSet.dataFourierPhase = nullptr; + } + if (dataSet.theory) { + delete dataSet.theory; + dataSet.theory = nullptr; + } + if (dataSet.theoryFourierRe) { + delete dataSet.theoryFourierRe; + dataSet.theoryFourierRe = nullptr; + } + if (dataSet.theoryFourierIm) { + delete dataSet.theoryFourierIm; + dataSet.theoryFourierIm = nullptr; + } + if (dataSet.theoryFourierPwr) { + delete dataSet.theoryFourierPwr; + dataSet.theoryFourierPwr = nullptr; + } + if (dataSet.theoryFourierPhase) { + delete dataSet.theoryFourierPhase; + dataSet.theoryFourierPhase = nullptr; + } + if (dataSet.diff) { + delete dataSet.diff; + dataSet.diff = nullptr; + } + if (dataSet.diffFourierRe) { + delete dataSet.diffFourierRe; + dataSet.diffFourierRe = nullptr; + } + if (dataSet.diffFourierIm) { + delete dataSet.diffFourierIm; + dataSet.diffFourierIm = nullptr; + } + if (dataSet.diffFourierPwr) { + delete dataSet.diffFourierPwr; + dataSet.diffFourierPwr = nullptr; + } + if (dataSet.diffFourierPhase) { + delete dataSet.diffFourierPhase; + dataSet.diffFourierPhase = nullptr; + } + if (dataSet.dataRange) { + delete dataSet.dataRange; + dataSet.dataRange = nullptr; + } +} + +//-------------------------------------------------------------------------- +// HandleDataSet (private) +//-------------------------------------------------------------------------- +/** + *

Generates the necessary histograms for plotting, starting from the pre-processed data. + * + * \param plotNo The number of the histo within the run list (fPlotNumber is the number of the plot BLOCK) + * \param runNo The number of the run + * \param data pre-processed data + */ +void PMusrViewModel::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) +{ + PMusrViewDataSet dataSet; + TH1F *dataHisto; + TH1F *theoHisto; + + TString name; + Double_t start; + Double_t end; + Double_t xmin, xmax, ymin, ymax; + Int_t size; + + InitDataSet(dataSet); + + // create plot range object for the data set + dataSet.dataRange = new PMusrViewPlotRange(); + + // dataHisto ------------------------------------------------------------- + // create histo specific infos + name = *fMsrHandler->GetMsrRunList()->at(runNo).GetRunName() + "_DataRunNo"; + name += static_cast(runNo); + name += "_"; + name += fPlotNumber; + start = data->GetDataTimeStart() - data->GetDataTimeStep()/2.0; + end = start + data->GetValue()->size()*data->GetDataTimeStep(); + size = data->GetValue()->size(); + dataSet.dataRange->SetXRange(start, end); // full possible range + // make sure that for asymmetry the y-range is initialized reasonably + if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) || (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR)) + dataSet.dataRange->SetYRange(-0.4, 0.4); + // extract necessary range information + if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() == 0) && + !fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { // no range information at all + if (fXRangePresent) { + if (fXmin > start) + fXmin = start; + if (fXmax < end) + fXmax = end; + } else { + fXRangePresent = true; + fXmin = start; + fXmax = end; + } + if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) || + (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR) || + (fMsrHandler->GetMsrRunList()->at(runNo).IsLifetimeCorrected())) { + fYRangePresent = true; + fYmin = -0.4; + fYmax = 0.4; + } + } + + // check if plot range is given in the msr-file, and if yes keep the values + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() == 1) { + // keep x-range + xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[0]; + xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[0]; + dataSet.dataRange->SetXRange(xmin, xmax); + // keep range information + fXRangePresent = true; + fXmin = xmin; + fXmax = xmax; + // check if y-range is given as well + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin.size() != 0) { + ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin[0]; + ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax[0]; + dataSet.dataRange->SetYRange(ymin, ymax); + // keep range information + fYRangePresent = true; + fYmin = ymin; + fYmax = ymax; + } + } + + // check if 'use_fit_ranges' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { + start = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(0); // needed to estimate size + if (start == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + start = fMsrHandler->GetMsrGlobal()->GetFitRange(0); + } + end = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(1); // needed to estimate size + if (end == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + end = fMsrHandler->GetMsrGlobal()->GetFitRange(1); + } + size = (Int_t) ((end - start) / data->GetDataTimeStep()) + 1; + start = data->GetDataTimeStart() + + (Int_t)((start - data->GetDataTimeStart())/data->GetDataTimeStep()) * data->GetDataTimeStep() - + data->GetDataTimeStep()/2.0; // closesd start value compatible with the user given + end = start + size * data->GetDataTimeStep(); // closesd end value compatible with the user given + dataSet.dataRange->SetXRange(start, end); + + // make sure that for asymmetry the y-range is initialized reasonably + if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) || + (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR) || + (fMsrHandler->GetMsrRunList()->at(runNo).IsLifetimeCorrected())) { + dataSet.dataRange->SetYRange(-0.4, 0.4); + } + + // keep range information + if (fXRangePresent) { + if (fXmin > start) + fXmin = start; + if (fXmax < end) + fXmax = end; + } else { + fXRangePresent = true; + fXmin = start; + fXmax = end; + } + // check if y-range is given as well + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin.size() != 0) { + ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin[0]; + ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax[0]; + dataSet.dataRange->SetYRange(ymin, ymax); + // keep range information + fYRangePresent = true; + fYmin = ymin; + fYmax = ymax; + } + } + + // check if 'sub_ranges' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() > 1) { + start = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo]; // needed to estimate size + end = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo]; // needed to estimate size + size = (Int_t) ((end - start) / data->GetDataTimeStep()) + 1; + start = data->GetDataTimeStart() + + (Int_t)((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep()) * data->GetDataTimeStep() - + data->GetDataTimeStep()/2.0; // closesd start value compatible with the user given + end = start + size * data->GetDataTimeStep(); // closesd end value compatible with the user given + dataSet.dataRange->SetXRange(start, end); + // keep range information + if (fXRangePresent) { + if (fXmin > start) + fXmin = start; + if (fXmax < end) + fXmax = end; + } else { + fXRangePresent = true; + fXmin = start; + fXmax = end; + } + + // check if y-range is given as well + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin.size() != 0) { + ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin[0]; + ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax[0]; + dataSet.dataRange->SetYRange(ymin, ymax); + // keep range information + fYRangePresent = true; + fYmin = ymin; + fYmax = ymax; + } + } + + // invoke histo + dataHisto = new TH1F(name, name, size, start, end); + + // fill histogram + // 1st calculate the bin-range according to the plot options + UInt_t startBin = 0; + UInt_t endBin = data->GetValue()->size(); + + // check if 'use_fit_range' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { + Double_t startFitRange = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(0); + if (startFitRange == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + startFitRange = fMsrHandler->GetMsrGlobal()->GetFitRange(0); + } + Double_t dval = (startFitRange - data->GetDataTimeStart())/data->GetDataTimeStep(); + if (dval < 0.0) { // make sure that startBin >= 0 + startBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin data < 0 for 'use_fit_range', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that startBin <= length of data vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',"; + std::cerr << std::endl << ">> will set it to data vector size" << std::endl << std::endl; + startBin = data->GetValue()->size(); + } else { + startBin = (UInt_t)dval; + } + + Double_t endFitRange = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(1); + if (endFitRange == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + endFitRange = fMsrHandler->GetMsrGlobal()->GetFitRange(1); + } + dval = (endFitRange - data->GetDataTimeStart())/data->GetDataTimeStep(); + if (dval < 0.0) { // make sure that endBin >= 0 + endBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin data < 0 for 'use_fit_range', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that endBin <= length of data vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',"; + std::cerr << std::endl << ">> will set it to data vector size" << std::endl << std::endl; + endBin = data->GetValue()->size(); + } else { + endBin = (UInt_t)dval; + } + } + + // check if 'sub_ranges' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() > 1) { + Double_t dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep(); + if (dval < 0.0) { // make sure that startBin >= 0 + startBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin data < 0 for 'sub_ranges', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that startBin <= length of data vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',"; + std::cerr << std::endl << ">> will set it to data vector size" << std::endl << std::endl; + startBin = data->GetValue()->size(); + } else { + startBin = (UInt_t)dval; + } + + dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep(); + if (dval < 0.0) { // make sure that endBin >= 0 + endBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin data < 0 for 'sub_ranges', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that endtBin <= length of data vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',"; + std::cerr << std::endl << ">> will set it to data vector size" << std::endl << std::endl; + endBin = data->GetValue()->size(); + } else { + endBin = (UInt_t)dval; + } + } + + for (UInt_t i=startBin; iSetBinContent(i-startBin+1, data->GetValue()->at(i)); + dataHisto->SetBinError(i-startBin+1, data->GetError()->at(i)); + } + + // set marker and line color + if (plotNo < fColorList.size()) { + dataHisto->SetMarkerColor(fColorList[plotNo]); + dataHisto->SetLineColor(fColorList[plotNo]); + } else { + TRandom rand(plotNo); + // the rgb values need to be drawn into separate variables first: the evaluation order of + // function arguments is unspecified in C++, hence calling rand.Integer() three times within + // the argument list would yield compiler dependent (r,g,b) permutations. + Int_t r = (Int_t)rand.Integer(255); + Int_t g = (Int_t)rand.Integer(255); + Int_t b = (Int_t)rand.Integer(255); + Int_t color = TColor::GetColor(r, g, b); + dataHisto->SetMarkerColor(color); + dataHisto->SetLineColor(color); + } + // set marker size + dataHisto->SetMarkerSize(1); + // set marker type + if (plotNo < fMarkerList.size()) { + dataHisto->SetMarkerStyle(fMarkerList[plotNo]); + } else { + TRandom rand(plotNo); + dataHisto->SetMarkerStyle(20+(Int_t)rand.Integer(10)); + } + + // theoHisto ------------------------------------------------------------- + // create histo specific infos + name = *fMsrHandler->GetMsrRunList()->at(runNo).GetRunName() + "_TheoRunNo"; + name += (Int_t)runNo; + name += "_"; + name += fPlotNumber; + start = data->GetTheoryTimeStart() - data->GetTheoryTimeStep()/2.0; + end = start + data->GetTheory()->size()*data->GetTheoryTimeStep(); + size = data->GetTheory()->size(); + + // check if 'use_fit_range' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { + start = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(0); // needed to estimate size + if (start == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + start = fMsrHandler->GetMsrGlobal()->GetFitRange(0); + } + end = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(1); // needed to estimate size + if (end == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + end = fMsrHandler->GetMsrGlobal()->GetFitRange(1); + } + size = (Int_t) ((end - start) / data->GetTheoryTimeStep()) + 1; + start = data->GetTheoryTimeStart() + + (Int_t)((start - data->GetTheoryTimeStart())/data->GetTheoryTimeStep()) * data->GetTheoryTimeStep() - + data->GetTheoryTimeStep()/2.0; // closesd start value compatible with the user given + end = start + size * data->GetTheoryTimeStep(); // closesd end value compatible with the user given + } + + // check if 'sub_ranges' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() > 1) { + start = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo]; // needed to estimate size + end = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo]; // needed to estimate size + size = (Int_t) ((end - start) / data->GetTheoryTimeStep()) + 1; + start = data->GetTheoryTimeStart() + + (Int_t)((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] - data->GetTheoryTimeStart())/data->GetTheoryTimeStep()) * data->GetTheoryTimeStep() - + data->GetTheoryTimeStep()/2.0; // closesd start value compatible with the user given + end = start + size * data->GetTheoryTimeStep(); // closesd end value compatible with the user given +} + + // invoke histo + theoHisto = new TH1F(name, name, size, start, end); + + // fill histogram + startBin = 0; + endBin = data->GetTheory()->size(); + + // check if 'use_fit_range' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { + Double_t startFitRange = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(0); + if (startFitRange == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + startFitRange = fMsrHandler->GetMsrGlobal()->GetFitRange(0); + } + Double_t dval = (startFitRange - data->GetDataTimeStart())/data->GetTheoryTimeStep(); + if (dval < 0.0) { // make sure that startBin >= 0 + startBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin theory < 0 for 'use_fit_range', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that startBin <= length of theory vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',"; + std::cerr << std::endl << ">> will set it to theory vector size" << std::endl << std::endl; + startBin = data->GetTheory()->size(); + } else { + startBin = (UInt_t)dval; + } + + Double_t endFitRange = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(1); + if (endFitRange == PMUSR_UNDEFINED) { // not given in the run block, try the global block entry + endFitRange = fMsrHandler->GetMsrGlobal()->GetFitRange(1); + } + dval = (endFitRange - data->GetDataTimeStart())/data->GetTheoryTimeStep(); + if (dval < 0.0) { // make sure that endBin >= 0 + endBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin theory < 0 for 'use_fit_range', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that endBin <= length of theory vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',"; + std::cerr << std::endl << ">> will set it to theory vector size" << std::endl << std::endl; + endBin = data->GetTheory()->size(); + } else { + endBin = (UInt_t)dval; + } + } + + // check if 'sub_ranges' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() > 1) { + startBin = (UInt_t)((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep()); + endBin = (UInt_t)((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep()); + + Double_t dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep(); + if (dval < 0.0) { // make sure that startBin >= 0 + startBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin theory < 0 for 'sub_ranges', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that startBin <= length of theory vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',"; + std::cerr << std::endl << ">> will set it to theory vector size" << std::endl << std::endl; + startBin = data->GetTheory()->size(); + } else { + startBin = (UInt_t)dval; + } + + dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep(); + if (dval < 0.0) { // make sure that endBin >= 0 + endBin = 0; + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin theory < 0 for 'sub_ranges', will set it to 0" << std::endl << std::endl; + } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that endtBin <= length of theory vector + std::cerr << std::endl << ">> PMusrViewModel::HandleDataSet(): **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',"; + std::cerr << std::endl << ">> will set it to theory vector size" << std::endl << std::endl; + endBin = data->GetTheory()->size(); + } else { + endBin = (UInt_t)dval; + } + } + + for (UInt_t i=startBin; iSetBinContent(i-startBin+1, data->GetTheory()->at(i)); + } + + // set the line color + if (plotNo < fColorList.size()) { + theoHisto->SetLineColor(fColorList[plotNo]); + } else { + TRandom rand(plotNo); + // the rgb values need to be drawn into separate variables first: the evaluation order of + // function arguments is unspecified in C++, hence calling rand.Integer() three times within + // the argument list would yield compiler dependent (r,g,b) permutations. + Int_t r = (Int_t)rand.Integer(255); + Int_t g = (Int_t)rand.Integer(255); + Int_t b = (Int_t)rand.Integer(255); + Int_t color = TColor::GetColor(r, g, b); + theoHisto->SetLineColor(color); + } + + // fill handler list ----------------------------------------------------- + dataSet.data = dataHisto; + dataSet.theory = theoHisto; + dataSet.diffFourierTag = 0; // not relevant at this point + + fData.push_back(dataSet); +} + +//-------------------------------------------------------------------------- +// HandleNonMusrDataSet (private) +//-------------------------------------------------------------------------- +/** + *

Generates the necessary error graphs for plotting, starting from the pre-processed data. + * + * \param plotNo The number of the histo within the run list (fPlotNumber is the number of the plot BLOCK) + * \param runNo The number of the run + * \param data pre-processed data + */ +void PMusrViewModel::HandleNonMusrDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) +{ + PMusrViewNonMusrDataSet dataSet; + TGraphErrors *dataHisto; + TGraphErrors *theoHisto; + + InitDataSet(dataSet); + + // create plot range object for the data set and fill it + dataSet.dataRange = new PMusrViewPlotRange(); + + // dataHisto ------------------------------------------------------------- + + // invoke graph + dataHisto = new TGraphErrors(data->GetX()->size()); + + // fill graph + for (UInt_t i=0; iGetX()->size(); i++) { + dataHisto->SetPoint(i, data->GetX()->at(i), data->GetValue()->at(i)); + dataHisto->SetPointError(i, 0.0, data->GetError()->at(i)); + } + + // set marker and line color + if (plotNo < fColorList.size()) { + dataHisto->SetMarkerColor(fColorList[plotNo]); + dataHisto->SetLineColor(fColorList[plotNo]); + } else { + TRandom rand(plotNo); + // the rgb values need to be drawn into separate variables first: the evaluation order of + // function arguments is unspecified in C++, hence calling rand.Integer() three times within + // the argument list would yield compiler dependent (r,g,b) permutations. + Int_t r = (Int_t)rand.Integer(255); + Int_t g = (Int_t)rand.Integer(255); + Int_t b = (Int_t)rand.Integer(255); + Int_t color = TColor::GetColor(r, g, b); + dataHisto->SetMarkerColor(color); + dataHisto->SetLineColor(color); + } + // set marker size + dataHisto->SetMarkerSize(1); + // set marker type + if (plotNo < fMarkerList.size()) { + dataHisto->SetMarkerStyle(fMarkerList[plotNo]); + } else { + TRandom rand(plotNo); + dataHisto->SetMarkerStyle(20+(Int_t)rand.Integer(10)); + } + + // theoHisto ------------------------------------------------------------- + + // invoke graph + theoHisto = new TGraphErrors(data->GetXTheory()->size()); + + // fill graph + for (UInt_t i=0; iGetXTheory()->size(); i++) { + theoHisto->SetPoint(i, data->GetXTheory()->at(i), data->GetTheory()->at(i)); + theoHisto->SetPointError(i, 0.0, 0.0); + } + + // set the line color + if (plotNo < fColorList.size()) { + theoHisto->SetLineColor(fColorList[plotNo]); + } else { + TRandom rand(plotNo); + // the rgb values need to be drawn into separate variables first: the evaluation order of + // function arguments is unspecified in C++, hence calling rand.Integer() three times within + // the argument list would yield compiler dependent (r,g,b) permutations. + Int_t r = (Int_t)rand.Integer(255); + Int_t g = (Int_t)rand.Integer(255); + Int_t b = (Int_t)rand.Integer(255); + Int_t color = TColor::GetColor(r, g, b); + theoHisto->SetLineColor(color); + } + + // fill handler list ----------------------------------------------------- + dataSet.data = dataHisto; + dataSet.theory = theoHisto; + dataSet.diffFourierTag = 0; // not relevant at this point + + // check the plot range options + Double_t xmin=0.0, xmax=0.0, ymin=0.0, ymax=0.0, x=0.0, y=0.0; + + // if no plot-range entry is present, initialize the plot range to the maximal possible given the data + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() == 0) { + dataSet.data->GetPoint(0, xmin, y); // get xmin + dataSet.data->GetPoint(dataSet.data->GetN()-1, xmax, y); // get xmax + dataSet.data->GetPoint(0, x, y); // init ymin/ymax + ymin = y; + ymax = y; + for (Int_t i=1; iGetN(); i++) { + dataSet.data->GetPoint(i, x, y); + if (y < ymin) + ymin = y; + if (y > ymax) + ymax = y; + } + Double_t dx = 0.025*(xmax-xmin); + Double_t dy = 0.025*(ymax-ymin); + dataSet.dataRange->SetXRange(xmin-dx, xmax+dx); + dataSet.dataRange->SetYRange(ymin-dy, ymax+dy); + } + + // check if plot range is given in the msr-file, and if yes keep the values + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() == 1) { + // keep x-range + xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[0]; + xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[0]; + dataSet.dataRange->SetXRange(xmin, xmax); + + // check if y-range is given as well + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin.size() != 0) { + ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin[0]; + ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax[0]; + dataSet.dataRange->SetYRange(ymin, ymax); + } + } + + // check if 'use_fit_range' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { + xmin = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(0); // needed to estimate size + xmax = fMsrHandler->GetMsrRunList()->at(runNo).GetFitRange(1); // needed to estimate size + dataSet.dataRange->SetXRange(xmin, xmax); + + // check if y-range is given as well + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin.size() != 0) { + ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin[0]; + ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax[0]; + dataSet.dataRange->SetYRange(ymin, ymax); + } + } + + // check if 'sub_ranges' plotting is whished + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() > 1) { + xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo]; // needed to estimate size + xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo]; // needed to estimate size + dataSet.dataRange->SetXRange(xmin, xmax); + + // check if y-range is given as well + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin.size() != 0) { + ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin[0]; + ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax[0]; + dataSet.dataRange->SetYRange(ymin, ymax); + } + } + + // keep maximal range of all plot present + fXRangePresent = true; + fYRangePresent = true; + if (plotNo == 0) { + fXmin = dataSet.dataRange->GetXmin(); + fXmax = dataSet.dataRange->GetXmax(); + fYmin = dataSet.dataRange->GetYmin(); + fYmax = dataSet.dataRange->GetYmax(); + } else { + if (fXmin > dataSet.dataRange->GetXmin()) + fXmin = dataSet.dataRange->GetXmin(); + if (fXmax < dataSet.dataRange->GetXmax()) + fXmax = dataSet.dataRange->GetXmax(); + if (fYmin > dataSet.dataRange->GetYmin()) + fYmin = dataSet.dataRange->GetYmin(); + if (fYmax < dataSet.dataRange->GetYmax()) + fYmax = dataSet.dataRange->GetYmax(); + } + + fNonMusrData.push_back(dataSet); +} + +//-------------------------------------------------------------------------- +// HandleDifference (private) +//-------------------------------------------------------------------------- +/** + *

Handles the calculation of the difference spectra (i.e. data-theory). + * It allocates the necessary objects if they are not already present. + */ +void PMusrViewModel::HandleDifference() +{ + // check if it is necessary to calculate diff data + if ((fPlotType != MSR_PLOT_NON_MUSR) && (fData[0].diff == nullptr)) { + TH1F *diffHisto; + TString name; + // loop over all histos + for (UInt_t i=0; iGetTitle()) + "_diff"; + diffHisto = new TH1F(name, name, fData[i].data->GetNbinsX(), + fData[i].data->GetXaxis()->GetXmin(), + fData[i].data->GetXaxis()->GetXmax()); + + // set marker and line color + diffHisto->SetMarkerColor(fData[i].data->GetMarkerColor()); + diffHisto->SetLineColor(fData[i].data->GetLineColor()); + // set marker size + diffHisto->SetMarkerSize(fData[i].data->GetMarkerSize()); + // set marker type + diffHisto->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + + // keep difference histo + fData[i].diff = diffHisto; + // calculate diff histo entry + double value; + for (Int_t j=1; jGetNbinsX()-1; j++) { + // set diff bin value + value = CalculateDiff(fData[i].data->GetBinCenter(j), + fData[i].data->GetBinContent(j), + fData[i].theory); + fData[i].diff->SetBinContent(j, value); + // set error diff bin value + value = fData[i].data->GetBinError(j); + fData[i].diff->SetBinError(j, value); + } + } + } else if ((fPlotType == MSR_PLOT_NON_MUSR) && (fNonMusrData[0].diff == nullptr)) { + TGraphErrors *diffHisto; + TString name; + // loop over all histos + for (UInt_t i=0; iGetN()); + + // create difference histos + name = TString(fNonMusrData[i].data->GetTitle()) + "_diff"; + diffHisto->SetNameTitle(name.Data(), name.Data()); + + // set marker and line color + diffHisto->SetMarkerColor(fNonMusrData[i].data->GetMarkerColor()); + diffHisto->SetLineColor(fNonMusrData[i].data->GetLineColor()); + // set marker size + diffHisto->SetMarkerSize(fNonMusrData[i].data->GetMarkerSize()); + // set marker type + diffHisto->SetMarkerStyle(fNonMusrData[i].data->GetMarkerStyle()); + + // keep difference histo + fNonMusrData[i].diff = diffHisto; + // calculate diff histo entry + double value; + double x, y; + for (Int_t j=0; jGetN(); j++) { + // set diff bin value + fNonMusrData[i].data->GetPoint(j, x, y); + value = CalculateDiff(x, y, fNonMusrData[i].theory); + fNonMusrData[i].diff->SetPoint(j, x, value); + // set error diff bin value + value = fNonMusrData[i].data->GetErrorY(j); + fNonMusrData[i].diff->SetPointError(j, 0.0, value); + } + } + } +} + +//-------------------------------------------------------------------------- +// HandleFourier (private) +//-------------------------------------------------------------------------- +/** + *

Handles the calculation of the Fourier transform. + * It allocates the necessary objects if they are not already present. + */ +void PMusrViewModel::HandleFourier(const PMusrViewMode &mode) +{ + Double_t re, im; + + // check if plot type is appropriate for fourier + if (fPlotType == MSR_PLOT_NON_MUSR) + return; + + // check if fourier needs to be calculated + if (fData[0].dataFourierRe == nullptr) { + Int_t bin; + double startTime = fXmin; + double endTime = fXmax; + if (fVisibleRangeSet && fHistoFrame) { // the user zoomed in, transform that window only + bin = fHistoFrame->GetXaxis()->GetFirst(); + startTime = fHistoFrame->GetBinLowEdge(bin); + bin = fHistoFrame->GetXaxis()->GetLast(); + endTime = fHistoFrame->GetBinLowEdge(bin)+fHistoFrame->GetBinWidth(bin); + } + // 1st: create the data/theory PFourier objects. FFTW plan creation is not + // thread-safe, hence this is done serially. + std::vector fourierData(fData.size(), nullptr); + std::vector fourierTheory(fData.size(), nullptr); + for (UInt_t i=0; iIsValid()) { + std::cerr << std::endl << ">> PMusrViewModel::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier data ..." << std::endl; + for (auto *ft : fourierData) delete ft; + for (auto *ft : fourierTheory) delete ft; + return; + } + + if (fTheoAsData) { // theory only at the data points + fourierTheory[i] = new PFourier(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower); + } else { + Int_t powerPad = fFourier.fFourierPower+5; // +5 means 8 times more points on theo (+3) + 4 times more points in fourier (+2) + fourierTheory[i] = new PFourier(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad); + } + if (!fourierTheory[i]->IsValid()) { + std::cerr << std::endl << ">> PMusrViewModel::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier theory ..." << std::endl; + for (auto *ft : fourierData) delete ft; + for (auto *ft : fourierTheory) delete ft; + return; + } + } + + // 2nd: execute the FFTs. Each data set is transformed independently, hence + // this can safely be parallelized. + Int_t noOfFT = static_cast(fData.size()); + #ifdef HAVE_GOMP + Int_t chunk = noOfFT/omp_get_num_procs(); + if (chunk < 1) + chunk = 1; + #pragma omp parallel for default(shared) schedule(dynamic,chunk) + #endif + for (Int_t i=0; iTransform(fFourier.fApodization); + fourierTheory[i]->Transform(fFourier.fApodization); + } + + // 3rd: extract the results into ROOT histograms and style them. ROOT + // object creation is not thread-safe, hence this is done serially. + for (UInt_t i=0; iGetBinWidth(1)/(endTime-startTime)); + // get real part of the data + fData[i].dataFourierRe = fourierData[i]->GetRealFourier(scale); + // get imaginary part of the data + fData[i].dataFourierIm = fourierData[i]->GetImaginaryFourier(scale); + // get power part of the data + fData[i].dataFourierPwr = fourierData[i]->GetPowerFourier(scale); + // get phase part of the data + fData[i].dataFourierPhase = fourierData[i]->GetPhaseFourier(); + + // set marker and line color + fData[i].dataFourierRe->SetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierRe->SetLineColor(fData[i].data->GetLineColor()); + fData[i].dataFourierIm->SetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierIm->SetLineColor(fData[i].data->GetLineColor()); + fData[i].dataFourierPwr->SetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPwr->SetLineColor(fData[i].data->GetLineColor()); + fData[i].dataFourierPhase->SetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPhase->SetLineColor(fData[i].data->GetLineColor()); + + // set marker size + fData[i].dataFourierRe->SetMarkerSize(1); + fData[i].dataFourierIm->SetMarkerSize(1); + fData[i].dataFourierPwr->SetMarkerSize(1); + fData[i].dataFourierPhase->SetMarkerSize(1); + + // set marker type + fData[i].dataFourierRe->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + fData[i].dataFourierIm->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + fData[i].dataFourierPwr->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + fData[i].dataFourierPhase->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + + scale = sqrt(fData[0].theory->GetBinWidth(1)/(endTime-startTime)*fData[0].theory->GetBinWidth(1)/fData[0].data->GetBinWidth(1)); + // get real part of the data + fData[i].theoryFourierRe = fourierTheory[i]->GetRealFourier(scale); + // get imaginary part of the data + fData[i].theoryFourierIm = fourierTheory[i]->GetImaginaryFourier(scale); + // get power part of the data + fData[i].theoryFourierPwr = fourierTheory[i]->GetPowerFourier(scale); + // get phase part of the data + fData[i].theoryFourierPhase = fourierTheory[i]->GetPhaseFourier(); + + // set line colors for the theory + fData[i].theoryFourierRe->SetLineColor(fData[i].theory->GetLineColor()); + fData[i].theoryFourierIm->SetLineColor(fData[i].theory->GetLineColor()); + fData[i].theoryFourierPwr->SetLineColor(fData[i].theory->GetLineColor()); + fData[i].theoryFourierPhase->SetLineColor(fData[i].theory->GetLineColor()); + + // clean up + delete fourierData[i]; + delete fourierTheory[i]; + } + + // phase opt. real FT requested initially in the msr-file, hence calculate it here + if (mode.GetFourierComponent() == EFourierComponent::kPhaseOptReal) { + CalcPhaseOptReFT(); + } + + // apply global phase if present + if (fFourier.fPhase.size() != 0.0) { + double cp; + double sp; + + fCurrentFourierPhase = fFourier.fPhase; + + for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set + // calculate new fourier data set value + re = fData[i].dataFourierRe->GetBinContent(j) * cp + fData[i].dataFourierIm->GetBinContent(j) * sp; + im = fData[i].dataFourierIm->GetBinContent(j) * cp - fData[i].dataFourierRe->GetBinContent(j) * sp; + // overwrite fourier data set value + fData[i].dataFourierRe->SetBinContent(j, re); + fData[i].dataFourierIm->SetBinContent(j, im); + } + } + if ((fData[i].theoryFourierRe != nullptr) && (fData[i].theoryFourierIm != nullptr)) { + for (Int_t j=0; jGetNbinsX(); j++) { // loop over a fourier data set + // calculate new fourier data set value + re = fData[i].theoryFourierRe->GetBinContent(j) * cp + fData[i].theoryFourierIm->GetBinContent(j) * sp; + im = fData[i].theoryFourierIm->GetBinContent(j) * cp - fData[i].theoryFourierRe->GetBinContent(j) * sp; + // overwrite fourier data set value + fData[i].theoryFourierRe->SetBinContent(j, re); + fData[i].theoryFourierIm->SetBinContent(j, im); + } + } + } + } + } +} + +//-------------------------------------------------------------------------- +// HandleDifferenceFourier (private) +//-------------------------------------------------------------------------- +/** + *

Handles the calculation of the Fourier transform of the difference spectra (i.e. data-theory). + * It allocates the necessary objects if they are not already present. + */ +void PMusrViewModel::HandleDifferenceFourier() +{ + // check if plot type is appropriate for fourier + if (fPlotType == MSR_PLOT_NON_MUSR) + return; + + // check if fourier needs to be calculated + if (fData[0].diffFourierRe == nullptr) { + // check if difference has been already calcualted, if not do it + if (fData[0].diff == nullptr) + HandleDifference(); + + // transform the range the user currently looks at + double startTime = fXmin; + double endTime = fXmax; + if (fVisibleRangeSet && fHistoFrame) { + Int_t bin = fHistoFrame->GetXaxis()->GetFirst(); + startTime = fHistoFrame->GetBinCenter(bin); + bin = fHistoFrame->GetXaxis()->GetLast(); + endTime = fHistoFrame->GetBinCenter(bin); + } + + // 1st: create the PFourier objects. FFTW plan creation is not + // thread-safe, hence this is done serially. + std::vector fourierData(fData.size(), nullptr); + for (UInt_t i=0; iIsValid()) { + std::cerr << std::endl << ">> PMusrViewModel::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier diff ..." << std::endl; + for (auto *ft : fourierData) delete ft; + return; + } + } + + // 2nd: execute the FFTs. Each data set is transformed independently, hence + // this can safely be parallelized. + Int_t noOfFT = static_cast(fData.size()); + #ifdef HAVE_GOMP + Int_t chunk = noOfFT/omp_get_num_procs(); + if (chunk < 1) + chunk = 1; + #pragma omp parallel for default(shared) schedule(dynamic,chunk) + #endif + for (Int_t i=0; iTransform(fFourier.fApodization); + } + + // 3rd: extract the results into ROOT histograms and style them. ROOT + // object creation is not thread-safe, hence this is done serially. + for (UInt_t i=0; iGetBinWidth(1)/(endTime-startTime)); + // get real part of the data + fData[i].diffFourierRe = fourierData[i]->GetRealFourier(scale); + // get imaginary part of the data + fData[i].diffFourierIm = fourierData[i]->GetImaginaryFourier(scale); + // get power part of the data + fData[i].diffFourierPwr = fourierData[i]->GetPowerFourier(scale); + // get phase part of the data + fData[i].diffFourierPhase = fourierData[i]->GetPhaseFourier(); + + // set marker and line color + fData[i].diffFourierRe->SetMarkerColor(fData[i].diff->GetMarkerColor()); + fData[i].diffFourierRe->SetLineColor(fData[i].diff->GetLineColor()); + fData[i].diffFourierIm->SetMarkerColor(fData[i].diff->GetMarkerColor()); + fData[i].diffFourierIm->SetLineColor(fData[i].diff->GetLineColor()); + fData[i].diffFourierPwr->SetMarkerColor(fData[i].diff->GetMarkerColor()); + fData[i].diffFourierPwr->SetLineColor(fData[i].diff->GetLineColor()); + fData[i].diffFourierPhase->SetMarkerColor(fData[i].diff->GetMarkerColor()); + fData[i].diffFourierPhase->SetLineColor(fData[i].diff->GetLineColor()); + + // set marker size + fData[i].diffFourierRe->SetMarkerSize(1); + fData[i].diffFourierIm->SetMarkerSize(1); + fData[i].diffFourierPwr->SetMarkerSize(1); + fData[i].diffFourierPhase->SetMarkerSize(1); + // set marker type + fData[i].diffFourierRe->SetMarkerStyle(fData[i].diff->GetMarkerStyle()); + fData[i].diffFourierIm->SetMarkerStyle(fData[i].diff->GetMarkerStyle()); + fData[i].diffFourierPwr->SetMarkerStyle(fData[i].diff->GetMarkerStyle()); + fData[i].diffFourierPhase->SetMarkerStyle(fData[i].diff->GetMarkerStyle()); + + // set diffFourierTag + fData[i].diffFourierTag = 1; // d-f + + // clean up + delete fourierData[i]; + } + + // apply phase + if (fFourier.fPhase.size() != 0.0) { + double re, im; + double cp; + double sp; + + fCurrentFourierPhase = fFourier.fPhase; + + for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set + // calculate new fourier data set value + re = fData[i].diffFourierRe->GetBinContent(j) * cp + fData[i].diffFourierIm->GetBinContent(j) * sp; + im = fData[i].diffFourierIm->GetBinContent(j) * cp - fData[i].diffFourierRe->GetBinContent(j) * sp; + // overwrite fourier data set value + fData[i].diffFourierRe->SetBinContent(j, re); + fData[i].diffFourierIm->SetBinContent(j, im); + } + } + } + } + } +} + +//-------------------------------------------------------------------------- +// HandleFourierDifference (private) +//-------------------------------------------------------------------------- +/** + *

Handles the calculation of the difference of the Fourier spectra. + * It allocates the necessary objects if they are not already present. + */ +void PMusrViewModel::HandleFourierDifference() +{ + // check if plot type is appropriate for fourier + if (fPlotType == MSR_PLOT_NON_MUSR) + return; + + // check if fourier needs to be calculated + if (fData[0].diffFourierRe == nullptr) { + // calculate all the Fourier differences + Double_t dval, dvalx; + TString name; + Int_t theoBin; + for (UInt_t i=0; iGetTitle()) + "_diff"; + fData[i].diffFourierRe = new TH1F(name, name, fData[i].dataFourierRe->GetNbinsX(), + fData[i].dataFourierRe->GetXaxis()->GetXmin(), + fData[i].dataFourierRe->GetXaxis()->GetXmax()); + + // imaginary part + name = TString(fData[i].dataFourierIm->GetTitle()) + "_diff"; + fData[i].diffFourierIm = new TH1F(name, name, fData[i].dataFourierIm->GetNbinsX(), + fData[i].dataFourierIm->GetXaxis()->GetXmin(), + fData[i].dataFourierIm->GetXaxis()->GetXmax()); + // power part + name = TString(fData[i].dataFourierPwr->GetTitle()) + "_diff"; + fData[i].diffFourierPwr = new TH1F(name, name, fData[i].dataFourierPwr->GetNbinsX(), + fData[i].dataFourierPwr->GetXaxis()->GetXmin(), + fData[i].dataFourierPwr->GetXaxis()->GetXmax()); + // phase part + name = TString(fData[i].dataFourierPhase->GetTitle()) + "_diff"; + fData[i].diffFourierPhase = new TH1F(name, name, fData[i].dataFourierPhase->GetNbinsX(), + fData[i].dataFourierPhase->GetXaxis()->GetXmin(), + fData[i].dataFourierPhase->GetXaxis()->GetXmax()); + + // phase optimized real part + if (fData[i].dataFourierPhaseOptReal != nullptr) { + name = TString(fData[i].dataFourierPhaseOptReal->GetTitle()) + "_diff"; + fData[i].diffFourierPhaseOptReal = new TH1F(name, name, fData[i].dataFourierPhaseOptReal->GetNbinsX(), + fData[i].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[i].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } + + // calculate difference + for (UInt_t j=1; jGetEntries(); j++) { + dvalx = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierRe->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierRe->GetBinContent(j) - fData[i].theoryFourierRe->GetBinContent(theoBin); + fData[i].diffFourierRe->SetBinContent(j, dval); + dvalx = fData[i].dataFourierIm->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierIm->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierIm->GetBinContent(j) - fData[i].theoryFourierIm->GetBinContent(theoBin); + fData[i].diffFourierIm->SetBinContent(j, dval); + dvalx = fData[i].dataFourierPwr->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierPwr->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierPwr->GetBinContent(j) - fData[i].theoryFourierPwr->GetBinContent(theoBin); + fData[i].diffFourierPwr->SetBinContent(j, dval); + dvalx = fData[i].dataFourierPhase->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierPhase->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierPhase->GetBinContent(j) - fData[i].theoryFourierPhase->GetBinContent(theoBin); + fData[i].diffFourierPhase->SetBinContent(j, dval); + if (fData[i].dataFourierPhaseOptReal != nullptr) { + dvalx = fData[i].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierPhaseOptReal->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierPhaseOptReal->GetBinContent(j) - fData[i].theoryFourierPhaseOptReal->GetBinContent(theoBin); + fData[i].diffFourierPhaseOptReal->SetBinContent(j, dval); + } + } + } + + for (UInt_t i=0; iSetMarkerColor(fData[i].dataFourierRe->GetMarkerColor()); + fData[i].diffFourierRe->SetLineColor(fData[i].dataFourierRe->GetLineColor()); + fData[i].diffFourierIm->SetMarkerColor(fData[i].dataFourierIm->GetMarkerColor()); + fData[i].diffFourierIm->SetLineColor(fData[i].dataFourierIm->GetLineColor()); + fData[i].diffFourierPwr->SetMarkerColor(fData[i].dataFourierPwr->GetMarkerColor()); + fData[i].diffFourierPwr->SetLineColor(fData[i].dataFourierPwr->GetLineColor()); + fData[i].diffFourierPhase->SetMarkerColor(fData[i].dataFourierPhase->GetMarkerColor()); + fData[i].diffFourierPhase->SetLineColor(fData[i].dataFourierPhase->GetLineColor()); + if (fData[i].dataFourierPhaseOptReal != nullptr) { + fData[i].diffFourierPhaseOptReal->SetMarkerColor(fData[i].dataFourierPhaseOptReal->GetMarkerColor()); + fData[i].diffFourierPhaseOptReal->SetLineColor(fData[i].dataFourierPhaseOptReal->GetLineColor()); + } + + // set marker size + fData[i].diffFourierRe->SetMarkerSize(1); + fData[i].diffFourierIm->SetMarkerSize(1); + fData[i].diffFourierPwr->SetMarkerSize(1); + fData[i].diffFourierPhase->SetMarkerSize(1); + if (fData[i].dataFourierPhaseOptReal != nullptr) { + fData[i].diffFourierPhaseOptReal->SetMarkerSize(1); + } + // set marker type + fData[i].diffFourierRe->SetMarkerStyle(fData[i].dataFourierRe->GetMarkerStyle()); + fData[i].diffFourierIm->SetMarkerStyle(fData[i].dataFourierIm->GetMarkerStyle()); + fData[i].diffFourierPwr->SetMarkerStyle(fData[i].dataFourierPwr->GetMarkerStyle()); + fData[i].diffFourierPhase->SetMarkerStyle(fData[i].dataFourierPhase->GetMarkerStyle()); + if (fData[i].dataFourierPhaseOptReal != nullptr) { + fData[i].diffFourierPhaseOptReal->SetMarkerStyle(fData[i].dataFourierPhaseOptReal->GetMarkerStyle()); + } + + // set diffFourierTag + fData[i].diffFourierTag = 2; // f-d + } + } +} + +//-------------------------------------------------------------------------- +// HandleAverage (private) +//-------------------------------------------------------------------------- +/** + *

Handles the calculation of the average of the ploted data. + * It allocates the necessary objects if they are not already present. + */ +void PMusrViewModel::HandleAverage() +{ + // check if plot type is appropriate for average + if (fPlotType == MSR_PLOT_NON_MUSR) + return; + + // in case there is still some average left over, cleanup first + if (fDataAvg.data != nullptr) { + CleanupAverage(); + } + + // create all the needed average data sets + TString name(""); + if (fData[0].data != nullptr) { + name = TString(fData[0].data->GetTitle()) + "_avg"; + fDataAvg.data = new TH1F(name, name, fData[0].data->GetNbinsX(), + fData[0].data->GetXaxis()->GetXmin(), + fData[0].data->GetXaxis()->GetXmax()); + } + if (fData[0].dataFourierRe != nullptr) { + name = TString(fData[0].dataFourierRe->GetTitle()) + "_avg"; + fDataAvg.dataFourierRe = new TH1F(name, name, fData[0].dataFourierRe->GetNbinsX(), + fData[0].dataFourierRe->GetXaxis()->GetXmin(), + fData[0].dataFourierRe->GetXaxis()->GetXmax()); + } + if (fData[0].dataFourierIm != nullptr) { + name = TString(fData[0].dataFourierIm->GetTitle()) + "_avg"; + fDataAvg.dataFourierIm = new TH1F(name, name, fData[0].dataFourierIm->GetNbinsX(), + fData[0].dataFourierIm->GetXaxis()->GetXmin(), + fData[0].dataFourierIm->GetXaxis()->GetXmax()); + } + if (fData[0].dataFourierPwr != nullptr) { + name = TString(fData[0].dataFourierPwr->GetTitle()) + "_avg"; + fDataAvg.dataFourierPwr = new TH1F(name, name, fData[0].dataFourierPwr->GetNbinsX(), + fData[0].dataFourierPwr->GetXaxis()->GetXmin(), + fData[0].dataFourierPwr->GetXaxis()->GetXmax()); + } + if (fData[0].dataFourierPhase != nullptr) { + name = TString(fData[0].dataFourierPhase->GetTitle()) + "_avg"; + fDataAvg.dataFourierPhase = new TH1F(name, name, fData[0].dataFourierPhase->GetNbinsX(), + fData[0].dataFourierPhase->GetXaxis()->GetXmin(), + fData[0].dataFourierPhase->GetXaxis()->GetXmax()); + } + if (fData[0].dataFourierPhaseOptReal != nullptr) { + name = TString(fData[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.dataFourierPhaseOptReal = new TH1F(name, name, fData[0].dataFourierPhaseOptReal->GetNbinsX(), + fData[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } + if (fData[0].theory != nullptr) { + name = TString(fData[0].theory->GetTitle()) + "_avg"; + fDataAvg.theory = new TH1F(name, name, fData[0].theory->GetNbinsX(), + fData[0].theory->GetXaxis()->GetXmin(), + fData[0].theory->GetXaxis()->GetXmax()); + } + if (fData[0].theoryFourierRe != nullptr) { + name = TString(fData[0].theoryFourierRe->GetTitle()) + "_avg"; + fDataAvg.theoryFourierRe = new TH1F(name, name, fData[0].theoryFourierRe->GetNbinsX(), + fData[0].theoryFourierRe->GetXaxis()->GetXmin(), + fData[0].theoryFourierRe->GetXaxis()->GetXmax()); + } + if (fData[0].theoryFourierIm != nullptr) { + name = TString(fData[0].theoryFourierIm->GetTitle()) + "_avg"; + fDataAvg.theoryFourierIm = new TH1F(name, name, fData[0].theoryFourierIm->GetNbinsX(), + fData[0].theoryFourierIm->GetXaxis()->GetXmin(), + fData[0].theoryFourierIm->GetXaxis()->GetXmax()); + } + if (fData[0].theoryFourierPwr != nullptr) { + name = TString(fData[0].theoryFourierPwr->GetTitle()) + "_avg"; + fDataAvg.theoryFourierPwr = new TH1F(name, name, fData[0].theoryFourierPwr->GetNbinsX(), + fData[0].theoryFourierPwr->GetXaxis()->GetXmin(), + fData[0].theoryFourierPwr->GetXaxis()->GetXmax()); + } + if (fData[0].theoryFourierPhase != nullptr) { + name = TString(fData[0].theoryFourierPhase->GetTitle()) + "_avg"; + fDataAvg.theoryFourierPhase = new TH1F(name, name, fData[0].theoryFourierPhase->GetNbinsX(), + fData[0].theoryFourierPhase->GetXaxis()->GetXmin(), + fData[0].theoryFourierPhase->GetXaxis()->GetXmax()); + } + if (fData[0].theoryFourierPhaseOptReal != nullptr) { + name = TString(fData[0].theoryFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.theoryFourierPhaseOptReal = new TH1F(name, name, fData[0].theoryFourierPhaseOptReal->GetNbinsX(), + fData[0].theoryFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].theoryFourierPhaseOptReal->GetXaxis()->GetXmax()); + } + if (fData[0].diff != nullptr) { + name = TString(fData[0].diff->GetTitle()) + "_avg"; + fDataAvg.diff = new TH1F(name, name, fData[0].diff->GetNbinsX(), + fData[0].diff->GetXaxis()->GetXmin(), + fData[0].diff->GetXaxis()->GetXmax()); + } + if (fData[0].diffFourierRe != nullptr) { + name = TString(fData[0].diffFourierRe->GetTitle()) + "_avg"; + fDataAvg.diffFourierRe = new TH1F(name, name, fData[0].diffFourierRe->GetNbinsX(), + fData[0].diffFourierRe->GetXaxis()->GetXmin(), + fData[0].diffFourierRe->GetXaxis()->GetXmax()); + } + if (fData[0].diffFourierIm != nullptr) { + name = TString(fData[0].diffFourierIm->GetTitle()) + "_avg"; + fDataAvg.diffFourierIm = new TH1F(name, name, fData[0].diffFourierIm->GetNbinsX(), + fData[0].diffFourierIm->GetXaxis()->GetXmin(), + fData[0].diffFourierIm->GetXaxis()->GetXmax()); + } + if (fData[0].diffFourierPwr != nullptr) { + name = TString(fData[0].diffFourierPwr->GetTitle()) + "_avg"; + fDataAvg.diffFourierPwr = new TH1F(name, name, fData[0].diffFourierPwr->GetNbinsX(), + fData[0].diffFourierPwr->GetXaxis()->GetXmin(), + fData[0].diffFourierPwr->GetXaxis()->GetXmax()); + } + if (fData[0].diffFourierPhase != nullptr) { + name = TString(fData[0].diffFourierPhase->GetTitle()) + "_avg"; + fDataAvg.diffFourierPhase = new TH1F(name, name, fData[0].diffFourierPhase->GetNbinsX(), + fData[0].diffFourierPhase->GetXaxis()->GetXmin(), + fData[0].diffFourierPhase->GetXaxis()->GetXmax()); + } + if (fData[0].diffFourierPhaseOptReal != nullptr) { + name = TString(fData[0].diffFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.diffFourierPhaseOptReal = new TH1F(name, name, fData[0].diffFourierPhaseOptReal->GetNbinsX(), + fData[0].diffFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].diffFourierPhaseOptReal->GetXaxis()->GetXmax()); + } + + // calculate all the average data sets + double dval; + if (fDataAvg.data != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.data->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.data->SetMarkerColor(fData[0].data->GetMarkerColor()); + fDataAvg.data->SetLineColor(fData[0].data->GetLineColor()); + fDataAvg.data->SetMarkerSize(fData[0].data->GetMarkerSize()); + fDataAvg.data->SetMarkerStyle(fData[0].data->GetMarkerStyle()); + } + if (fDataAvg.dataFourierRe != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierRe->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierRe->SetMarkerColor(fData[0].dataFourierRe->GetMarkerColor()); + fDataAvg.dataFourierRe->SetLineColor(fData[0].dataFourierRe->GetLineColor()); + fDataAvg.dataFourierRe->SetMarkerSize(fData[0].dataFourierRe->GetMarkerSize()); + fDataAvg.dataFourierRe->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle()); + } + if (fDataAvg.dataFourierIm != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierIm->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierIm->SetMarkerColor(fData[0].dataFourierIm->GetMarkerColor()); + fDataAvg.dataFourierIm->SetLineColor(fData[0].dataFourierIm->GetLineColor()); + fDataAvg.dataFourierIm->SetMarkerSize(fData[0].dataFourierIm->GetMarkerSize()); + fDataAvg.dataFourierIm->SetMarkerStyle(fData[0].dataFourierIm->GetMarkerStyle()); + } + if (fDataAvg.dataFourierPwr != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierPwr->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierPwr->SetMarkerColor(fData[0].dataFourierPwr->GetMarkerColor()); + fDataAvg.dataFourierPwr->SetLineColor(fData[0].dataFourierPwr->GetLineColor()); + fDataAvg.dataFourierPwr->SetMarkerSize(fData[0].dataFourierPwr->GetMarkerSize()); + fDataAvg.dataFourierPwr->SetMarkerStyle(fData[0].dataFourierPwr->GetMarkerStyle()); + } + if (fDataAvg.dataFourierPhase != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierPhase->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierPhase->SetMarkerColor(fData[0].dataFourierPhase->GetMarkerColor()); + fDataAvg.dataFourierPhase->SetLineColor(fData[0].dataFourierPhase->GetLineColor()); + fDataAvg.dataFourierPhase->SetMarkerSize(fData[0].dataFourierPhase->GetMarkerSize()); + fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); + } + if (fDataAvg.dataFourierPhaseOptReal != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierPhaseOptReal->SetMarkerColor(fData[0].dataFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.dataFourierPhaseOptReal->SetLineColor(fData[0].dataFourierPhaseOptReal->GetLineColor()); + fDataAvg.dataFourierPhaseOptReal->SetMarkerSize(fData[0].dataFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.dataFourierPhaseOptReal->SetMarkerStyle(fData[0].dataFourierPhaseOptReal->GetMarkerStyle()); + } + if (fDataAvg.theory != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theory->SetBinContent(i, dval/fData.size()); + } + fDataAvg.theory->SetLineColor(fData[0].theory->GetLineColor()); + } + if (fDataAvg.theoryFourierRe != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierRe->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierRe->SetMarkerColor(fData[0].theoryFourierRe->GetMarkerColor()); + fDataAvg.theoryFourierRe->SetLineColor(fData[0].theoryFourierRe->GetLineColor()); + fDataAvg.theoryFourierRe->SetMarkerSize(fData[0].theoryFourierRe->GetMarkerSize()); + fDataAvg.theoryFourierRe->SetMarkerStyle(fData[0].theoryFourierRe->GetMarkerStyle()); + } + if (fDataAvg.theoryFourierIm != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierIm->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierIm->SetMarkerColor(fData[0].theoryFourierIm->GetMarkerColor()); + fDataAvg.theoryFourierIm->SetLineColor(fData[0].theoryFourierIm->GetLineColor()); + fDataAvg.theoryFourierIm->SetMarkerSize(fData[0].theoryFourierIm->GetMarkerSize()); + fDataAvg.theoryFourierIm->SetMarkerStyle(fData[0].theoryFourierIm->GetMarkerStyle()); + } + if (fDataAvg.theoryFourierPwr != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierPwr->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierPwr->SetMarkerColor(fData[0].theoryFourierPwr->GetMarkerColor()); + fDataAvg.theoryFourierPwr->SetLineColor(fData[0].theoryFourierPwr->GetLineColor()); + fDataAvg.theoryFourierPwr->SetMarkerSize(fData[0].theoryFourierPwr->GetMarkerSize()); + fDataAvg.theoryFourierPwr->SetMarkerStyle(fData[0].theoryFourierPwr->GetMarkerStyle()); + } + if (fDataAvg.theoryFourierPhase != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierPhase->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierPhase->SetMarkerColor(fData[0].theoryFourierPhase->GetMarkerColor()); + fDataAvg.theoryFourierPhase->SetLineColor(fData[0].theoryFourierPhase->GetLineColor()); + fDataAvg.theoryFourierPhase->SetMarkerSize(fData[0].theoryFourierPhase->GetMarkerSize()); + fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle()); + } + if (fDataAvg.theoryFourierPhaseOptReal != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierPhaseOptReal->SetMarkerColor(fData[0].theoryFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.theoryFourierPhaseOptReal->SetLineColor(fData[0].theoryFourierPhaseOptReal->GetLineColor()); + fDataAvg.theoryFourierPhaseOptReal->SetMarkerSize(fData[0].theoryFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.theoryFourierPhaseOptReal->SetMarkerStyle(fData[0].theoryFourierPhaseOptReal->GetMarkerStyle()); + } + if (fDataAvg.diff != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.diff->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.diff->SetMarkerColor(fData[0].diff->GetMarkerColor()); + fDataAvg.diff->SetLineColor(fData[0].diff->GetLineColor()); + fDataAvg.diff->SetMarkerSize(fData[0].diff->GetMarkerSize()); + fDataAvg.diff->SetMarkerStyle(fData[0].diff->GetMarkerStyle()); + } + if (fDataAvg.diffFourierRe != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.diffFourierRe->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierRe->SetMarkerColor(fData[0].diffFourierRe->GetMarkerColor()); + fDataAvg.diffFourierRe->SetLineColor(fData[0].diffFourierRe->GetLineColor()); + fDataAvg.diffFourierRe->SetMarkerSize(fData[0].diffFourierRe->GetMarkerSize()); + fDataAvg.diffFourierRe->SetMarkerStyle(fData[0].diffFourierRe->GetMarkerStyle()); + } + if (fDataAvg.diffFourierIm != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.diffFourierIm->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierIm->SetMarkerColor(fData[0].diffFourierIm->GetMarkerColor()); + fDataAvg.diffFourierIm->SetLineColor(fData[0].diffFourierIm->GetLineColor()); + fDataAvg.diffFourierIm->SetMarkerSize(fData[0].diffFourierIm->GetMarkerSize()); + fDataAvg.diffFourierIm->SetMarkerStyle(fData[0].diffFourierIm->GetMarkerStyle()); + } + if (fDataAvg.diffFourierPwr != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.diffFourierPwr->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierPwr->SetMarkerColor(fData[0].diffFourierPwr->GetMarkerColor()); + fDataAvg.diffFourierPwr->SetLineColor(fData[0].diffFourierPwr->GetLineColor()); + fDataAvg.diffFourierPwr->SetMarkerSize(fData[0].diffFourierPwr->GetMarkerSize()); + fDataAvg.diffFourierPwr->SetMarkerStyle(fData[0].diffFourierPwr->GetMarkerStyle()); + } + if (fDataAvg.diffFourierPhase != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierPhase->SetMarkerColor(fData[0].diffFourierPhase->GetMarkerColor()); + fDataAvg.diffFourierPhase->SetLineColor(fData[0].diffFourierPhase->GetLineColor()); + fDataAvg.diffFourierPhase->SetMarkerSize(fData[0].diffFourierPhase->GetMarkerSize()); + fDataAvg.diffFourierPhase->SetMarkerStyle(fData[0].diffFourierPhase->GetMarkerStyle()); + } + if (fDataAvg.diffFourierPhaseOptReal != nullptr) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.diffFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierPhaseOptReal->SetMarkerColor(fData[0].dataFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.diffFourierPhaseOptReal->SetLineColor(fData[0].dataFourierPhaseOptReal->GetLineColor()); + fDataAvg.diffFourierPhaseOptReal->SetMarkerSize(fData[0].dataFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.diffFourierPhaseOptReal->SetMarkerStyle(fData[0].dataFourierPhaseOptReal->GetMarkerStyle()); + } +} + +//-------------------------------------------------------------------------- +// CleanupDifference (private) +//-------------------------------------------------------------------------- +/** + *

Cleans up (deallocate) difference data. + */ +void PMusrViewModel::CleanupDifference() +{ + for (UInt_t i=0; iCleans up (deallocate) Fourier transform data. + */ +void PMusrViewModel::CleanupFourier() +{ + for (UInt_t i=0; iCleans up (deallocate) Fourier difference spectra. + */ +void PMusrViewModel::CleanupFourierDifference() +{ + for (UInt_t i=0; iCleans up (deallocate) averaged data set. + */ +void PMusrViewModel::CleanupAverage() +{ + if (fDataAvg.data != nullptr) { + delete fDataAvg.data; + fDataAvg.data = nullptr; + } + if (fDataAvg.dataFourierRe != nullptr) { + delete fDataAvg.dataFourierRe; + fDataAvg.dataFourierRe = nullptr; + } + if (fDataAvg.dataFourierIm != nullptr) { + delete fDataAvg.dataFourierIm; + fDataAvg.dataFourierIm = nullptr; + } + if (fDataAvg.dataFourierPwr != nullptr) { + delete fDataAvg.dataFourierPwr; + fDataAvg.dataFourierPwr = nullptr; + } + if (fDataAvg.dataFourierPhase != nullptr) { + delete fDataAvg.dataFourierPhase; + fDataAvg.dataFourierPhase = nullptr; + } + if (fDataAvg.dataFourierPhaseOptReal != nullptr) { + delete fDataAvg.dataFourierPhaseOptReal; + fDataAvg.dataFourierPhaseOptReal = nullptr; + } + if (fDataAvg.theory != nullptr) { + delete fDataAvg.theory; + fDataAvg.theory = nullptr; + } + if (fDataAvg.theoryFourierRe != nullptr) { + delete fDataAvg.theoryFourierRe; + fDataAvg.theoryFourierRe = nullptr; + } + if (fDataAvg.theoryFourierIm != nullptr) { + delete fDataAvg.theoryFourierIm; + fDataAvg.theoryFourierIm = nullptr; + } + if (fDataAvg.theoryFourierPwr != nullptr) { + delete fDataAvg.theoryFourierPwr; + fDataAvg.theoryFourierPwr = nullptr; + } + if (fDataAvg.theoryFourierPhase != nullptr) { + delete fDataAvg.theoryFourierPhase; + fDataAvg.theoryFourierPhase = nullptr; + } + if (fDataAvg.theoryFourierPhaseOptReal != nullptr) { + delete fDataAvg.theoryFourierPhaseOptReal; + fDataAvg.theoryFourierPhaseOptReal = nullptr; + } + if (fDataAvg.diff != nullptr) { + delete fDataAvg.diff; + fDataAvg.diff = nullptr; + } + if (fDataAvg.diffFourierRe != nullptr) { + delete fDataAvg.diffFourierRe; + fDataAvg.diffFourierRe = nullptr; + } + if (fDataAvg.diffFourierIm != nullptr) { + delete fDataAvg.diffFourierIm; + fDataAvg.diffFourierIm = nullptr; + } + if (fDataAvg.diffFourierPwr != nullptr) { + delete fDataAvg.diffFourierPwr; + fDataAvg.diffFourierPwr = nullptr; + } + if (fDataAvg.diffFourierPhase != nullptr) { + delete fDataAvg.diffFourierPhase; + fDataAvg.diffFourierPhase = nullptr; + } + if (fDataAvg.diffFourierPhaseOptReal != nullptr) { + delete fDataAvg.diffFourierPhaseOptReal; + fDataAvg.diffFourierPhaseOptReal = nullptr; + } +} + +//-------------------------------------------------------------------------- +// CalculateDiff (private) +//-------------------------------------------------------------------------- +/** + * @brief PMusrViewModel::CalcPhaseOptReFT + */ +void PMusrViewModel::CalcPhaseOptReFT() +{ + Double_t min = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; + Double_t max = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; + + if ((min == -1.0) && (max == -1.0)) { + if ((fFourier.fPlotRange[0] != -1) && (fFourier.fPlotRange[1] != -1)) { + min = fFourier.fPlotRange[0]; + max = fFourier.fPlotRange[1]; + } else { + min = fData[0].dataFourierRe->GetBinLowEdge(1); + max = fData[0].dataFourierRe->GetBinLowEdge(fData[0].dataFourierRe->GetNbinsX())+fData[0].dataFourierRe->GetBinWidth(1); + } + } + + PDoubleVector phaseParam; + Char_t hName[1024]; + Double_t ph, re; + + for (UInt_t i=0; iSetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPhaseOptReal->SetLineColor(fData[i].data->GetLineColor()); + // set marker size + fData[i].dataFourierPhaseOptReal->SetMarkerSize(1); + // set marker type + fData[i].dataFourierPhaseOptReal->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + + // handle Fourier theory part + // clone theory Re FT + strcpy(hName, fData[i].theoryFourierPhase->GetName()); + strcat(hName, "_Opt_Real"); + fData[i].theoryFourierPhaseOptReal = (TH1F*) fData[i].theoryFourierRe->Clone(hName); + + // rotate the theory according to the optimized phase parameters + // first find minBin for min of the phase correction + Int_t minBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindFixBin(min); + Int_t maxBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindFixBin(max); + + for (Int_t j=1; jGetNbinsX(); j++) { + ph = phaseParam[0] + phaseParam[1] * (Double_t)(j-minBin+1) / (Double_t)(maxBin-minBin); + re = fData[i].theoryFourierRe->GetBinContent(j) * cos(ph) - fData[i].theoryFourierIm->GetBinContent(j) * sin(ph); + fData[i].theoryFourierPhaseOptReal->SetBinContent(j, re); + } + // set line colors for the theory + fData[i].theoryFourierPhaseOptReal->SetLineColor(fData[i].theory->GetLineColor()); + } +} + +//-------------------------------------------------------------------------- +// CalculateDiff (private) +//-------------------------------------------------------------------------- +/** + *

Calculates the difference between data and theory for histograms. + * + * return: + * - (data - theory) value + * + * \param x x-value of the data + * \param y y-value of the data + * \param theo theory histogram + */ +Double_t PMusrViewModel::CalculateDiff(const Double_t x, const Double_t y, TH1F *theo) +{ + Int_t bin = theo->FindBin(x); + + return y - theo->GetBinContent(bin); +} + +//-------------------------------------------------------------------------- +// CalculateDiff (private) +//-------------------------------------------------------------------------- +/** + *

Calculates the difference between data and theory for error graphs. + * + * return: + * - (data - theory) value + * + * \param x x-value of the data + * \param y y-value of the data + * \param theo theory error graphs + */ +Double_t PMusrViewModel::CalculateDiff(const Double_t x, const Double_t y, TGraphErrors *theo) +{ + Int_t bin = 0; + Double_t xVal, yVal; + + bin = FindBin(x, theo); + + theo->GetPoint(bin, xVal, yVal); + + return y - yVal; +} + +//-------------------------------------------------------------------------- +// FindBin (private) +//-------------------------------------------------------------------------- +/** + *

Analog to FindBin for histograms (TH1F) but here for TGraphErrors. + * + * return: + * - bin closest to a given x value. + * + * \param x x-value of the data + * \param graph TGraphErrors which should be searched + */ +Int_t PMusrViewModel::FindBin(const Double_t x, TGraphErrors *graph) +{ + Int_t i, bin = 0; + Double_t *xTheo = graph->GetX(); + + // find proper bin of the graph + for (i=0; iGetN(); i++) { + if (*(xTheo+i) >= x) { + bin = i; + break; + } + } + // in case it is the last point + if (i == graph->GetN()) { + bin = i; + } + + return bin; +} + +//-------------------------------------------------------------------------- +// GetMaximum (private) +//-------------------------------------------------------------------------- +/** + *

returns the maximum of a histogram in the range [xmin, xmax]. + * If xmin = xmax = -1.0, the whole histogram range is used. + * + * return: + * - maximum, or 0.0 if the histo pointer is the null pointer. + * + * \param histo pointer of the histogram + * \param xmin lower edge for the search interval. + * \param xmax upper edge for the search interval. + */ +Double_t PMusrViewModel::GetMaximum(TH1F* histo, Double_t xmin, Double_t xmax) +{ + if (histo == nullptr) + return 0.0; + + Int_t start=0, end=0; + if (xmin == xmax) { + start = 1; + end = histo->GetNbinsX(); + } else { + start = histo->FindBin(xmin); + if ((start==0) || (start==histo->GetNbinsX()+1)) // underflow/overflow + start = 1; + end = histo->FindBin(xmax); + if ((end==0) || (end==histo->GetNbinsX()+1)) // underflow/overflow + end = histo->GetNbinsX(); + } + + Double_t max = histo->GetBinContent(start); + Double_t binContent; + for (Int_t i=start; iGetBinContent(i); + if (max < binContent) + max = binContent; + } + + return max; +} + +//-------------------------------------------------------------------------- +// GetMinimum (private) +//-------------------------------------------------------------------------- +/** + *

returns the minimum of a histogram in the range [xmin, xmax]. + * If xmin = xmax = -1.0, the whole histogram range is used. + * + * return: + * - minimum, or 0.0 if the histo pointer is the null pointer. + * + * \param histo pointer of the histogram + * \param xmin lower edge for the search interval. + * \param xmax upper edge for the search interval. + */ +Double_t PMusrViewModel::GetMinimum(TH1F* histo, Double_t xmin, Double_t xmax) +{ + if (histo == nullptr) + return 0.0; + + Int_t start=0, end=0; + if (xmin == xmax) { + start = 1; + end = histo->GetNbinsX(); + } else { + start = histo->FindBin(xmin); + if ((start==0) || (start==histo->GetNbinsX()+1)) // underflow/overflow + start = 1; + end = histo->FindBin(xmax); + if ((end==0) || (end==histo->GetNbinsX()+1)) // underflow/overflow + end = histo->GetNbinsX(); + } + + Double_t min = histo->GetBinContent(start); + Double_t binContent; + for (Int_t i=start; iGetBinContent(i); + if (min > binContent) + min = binContent; + } + + return min; +} + +//-------------------------------------------------------------------------- +// GetMaximum (private) +//-------------------------------------------------------------------------- +/** + *

returns the maximum of a TGraphErrors object in the range [xmin, xmax]. + * If xmin = xmax = -1.0, the whole histogram range is used. + * + * return: + * - maximum, or 0.0 if the histo pointer is the null pointer. + * + * \param graph pointer of the histogram + * \param xmin lower edge for the search interval. + * \param xmax upper edge for the search interval. + */ +Double_t PMusrViewModel::GetMaximum(TGraphErrors* graph, Double_t xmin, Double_t xmax) +{ + if (graph == nullptr) + return 0.0; + + Double_t x, y; + if (xmin == xmax) { + graph->GetPoint(0, x, y); + xmin = x; + graph->GetPoint(graph->GetN()-1, x, y); + xmax = x; + } + + graph->GetPoint(0, x, y); + Double_t max = y; + for (Int_t i=0; iGetN(); i++) { + graph->GetPoint(i, x, y); + if ((x >= xmin) && (x <= xmax)) { + if (y > max) + max = y; + } + } + + return max; +} + +//-------------------------------------------------------------------------- +// GetMinimum (private) +//-------------------------------------------------------------------------- +/** + *

returns the minimum of a TGraphErrors object in the range [xmin, xmax]. + * If xmin = xmax = -1.0, the whole histogram range is used. + * + * return: + * - minimum, or 0.0 if the histo pointer is the null pointer. + * + * \param graph pointer of the histogram + * \param xmin lower edge for the search interval. + * \param xmax upper edge for the search interval. + */ +Double_t PMusrViewModel::GetMinimum(TGraphErrors* graph, Double_t xmin, Double_t xmax) +{ + if (graph == nullptr) + return 0.0; + + Double_t x, y; + if (xmin == xmax) { + graph->GetPoint(0, x, y); + xmin = x; + graph->GetPoint(graph->GetN()-1, x, y); + xmax = x; + } + + graph->GetPoint(0, x, y); + Double_t min = y; + for (Int_t i=0; iGetN(); i++) { + graph->GetPoint(i, x, y); + if ((x >= xmin) && (x <= xmax)) { + if (y < min) + min = y; + } + } + + return min; +} + +//-------------------------------------------------------------------------- +// IsScaleN0AndBkg (private) +//-------------------------------------------------------------------------- +/** + *

Checks if N0/Bkg normalization to 1/ns is whished. The default is yes, since most of the users want to have it that way. + * To overwrite this, one should add the line 'SCALE_N0_BKG FALSE' to the command block of the msr-file. + * + * return: + * - true, if scaling of N0 and Bkg to 1/ns is whished + * - false, otherwise + * + * \param histoNo forward histogram number of the run + */ +Bool_t PMusrViewModel::IsScaleN0AndBkg() +{ + Bool_t willScale = true; + + PMsrLines *cmd = fMsrHandler->GetMsrCommands(); + for (UInt_t i=0; isize(); i++) { + if (cmd->at(i).fLine.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) { + std::vector tokens = PStringUtils::Split(cmd->at(i).fLine.Data(), " \t"); + if (tokens.size() != 2) { + std::cerr << std::endl << ">> PRunSingleHisto::IsScaleN0AndBkg(): **WARNING** Found uncorrect 'SCALE_N0_BKG' command, will ignore it."; + std::cerr << std::endl << ">> Allowed commands: SCALE_N0_BKG TRUE | FALSE" << std::endl; + return willScale; + } + if (PStringUtils::IsEqualNoCase(tokens[1], "FALSE")) { + willScale = false; + } + } + } + + return willScale; +} + +//-------------------------------------------------------------------------- +// GetNeededAccuracy (private) +//-------------------------------------------------------------------------- +/** + *

Calculates the needed accuracy of the parameter value based on the given errors. + * + * return: + * - needed accuracy + * + * \param param fit parameter with its additional informations, like errors etc. + */ +UInt_t PMusrViewModel::GetNeededAccuracy(PMsrParamStructure param) +{ + const UInt_t precLimit = 6; + UInt_t decimalPoint = 0; + UInt_t accuracy = 6; + + if (param.fStep == 0.0) { // check if fit parameter is a constant, i.e. step==0 + char str[128]; + + snprintf(str, sizeof(str), "%lf", param.fValue); + + // find decimal point + for (UInt_t i=0; i=0; i--) { + if (str[i] != '0') { + if (((UInt_t)i-decimalPoint) < precLimit) + accuracy = (UInt_t)i-decimalPoint; + else + accuracy = precLimit; + break; + } + } + + } else if ((param.fStep != 0) && !param.fPosErrorPresent) { // check if no positive error is given step!=0 but fPosErrorPresent==false + + for (UInt_t i=0; i accuracy) + accuracy = accuracy2; + + if (accuracy+1 <= precLimit) + accuracy += 1; + } + + return accuracy; +} + +//-------------------------------------------------------------------------- +// GetInterpolatedValue (private) +//-------------------------------------------------------------------------- +/** + *

search for xVal in histo. If xVal is not found exactly, interpolate and + * return the interpolated y-value. + * + * return: + * - interpolated value if xVal is within histo range, 0 otherwise. + * + * \param histo pointer of the histogram + * \param xVal x-value to be looked for + */ +Double_t PMusrViewModel::GetInterpolatedValue(TH1F* histo, Double_t xVal) +{ + if (histo == nullptr) + return 0.0; + + Int_t idx = histo->FindBin(xVal); + + // make sure idx is within range + if ((idx < 1) || (idx > histo->GetNbinsX())) + return 0.0; + + // make sure the lower bound idx is found. This is needed since + // FindBin rounds towards the closer idx. + if (histo->GetBinCenter(idx) > xVal) + --idx; + + Double_t x0, x1, y0, y1; + x0 = histo->GetBinCenter(idx); + x1 = histo->GetBinCenter(idx+1); + y0 = histo->GetBinContent(idx); + y1 = histo->GetBinContent(idx+1); + + return (y1-y0)*(xVal-x0)/(x1-x0)+y0; +} + +//-------------------------------------------------------------------------- +// GetExportDataSet (private) +//-------------------------------------------------------------------------- +/** + *

extract data for export. + * + * \param data + * \param xmin + * \param xmax + * \param dumpData + * \param hasError + */ +void PMusrViewModel::GetExportDataSet(const TH1F *data, const Double_t xmin, const Double_t xmax, + PMusrViewAsciiDumpVector &dumpData, const Bool_t hasError) +{ + PMusrViewAsciiDump dump; + Double_t x=0.0; + + // go through all difference data bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get time/freq + x = data->GetBinCenter(j); + // check if x is in the current range + if ((x >= xmin) && (x <= xmax)) { + dump.dataX.push_back(x); + dump.data.push_back(data->GetBinContent(j)); + if (hasError) + dump.dataErr.push_back(data->GetBinError(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpData.push_back(dump); +} + +} // namespace musrview diff --git a/src/classes/PMusrWebCanvas.cpp b/src/classes/PMusrWebCanvas.cpp new file mode 100644 index 00000000..3c772f34 --- /dev/null +++ b/src/classes/PMusrWebCanvas.cpp @@ -0,0 +1,1117 @@ +/*************************************************************************** + + PMusrWebCanvas.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2026 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 + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "PMusrWebCanvas.h" + +namespace musrview { + +//-------------------------------------------------------------------------- +// local helpers +//-------------------------------------------------------------------------- +namespace { + +//------------------------------------------------------------------------ +/** + *

Locates the directory holding the musrview web assets (html/css/js). + * + *

The search order is: the environment variable MUSRFIT_WEB_DIR (used during + * development), then the directory configured at build time, then a couple of + * standard locations. + * + * \return the directory, or an empty string if none was found + */ +//------------------------------------------------------------------------ +/** + *

Converts a ROOT colour index into an HTML colour. + * + * \param colorIndex ROOT colour index + * + * \return the colour as "#rrggbb" + */ +std::string RootColorToHtml(Int_t colorIndex) +{ + TColor *col = gROOT->GetColor(colorIndex); + if (col == nullptr) + return "#000000"; + + Char_t buf[16]; + snprintf(buf, sizeof(buf), "#%02x%02x%02x", + static_cast(255*col->GetRed() + 0.5), + static_cast(255*col->GetGreen() + 0.5), + static_cast(255*col->GetBlue() + 0.5)); + + return std::string(buf); +} + +//------------------------------------------------------------------------ +/** + *

Maps a ROOT marker style onto the shape the web front end can draw. + * + * \param markerStyle ROOT marker style + * \param shape where the shape name is returned + * \param filled where the fill flag is returned + */ +void RootMarkerToShape(Int_t markerStyle, std::string &shape, Bool_t &filled) +{ + filled = kFALSE; + + switch (markerStyle) { + case 1: case 6: case 7: case 8: shape = "circle"; filled = kTRUE; break; + case 2: case 28: case 34: shape = "cross"; filled = (markerStyle == 34); break; + case 3: case 31: shape = "star"; break; + case 4: case 24: shape = "circle"; break; + case 5: shape = "xcross"; break; + case 20: shape = "circle"; filled = kTRUE; break; + case 21: shape = "square"; filled = kTRUE; break; + case 22: shape = "triangle-up"; filled = kTRUE; break; + case 23: shape = "triangle-down"; filled = kTRUE; break; + case 25: shape = "square"; break; + case 26: shape = "triangle-up"; break; + case 27: shape = "diamond"; break; + case 29: shape = "star"; filled = kTRUE; break; + case 30: shape = "star"; break; + case 32: shape = "triangle-down"; break; + case 33: shape = "diamond"; filled = kTRUE; break; + default: shape = "circle"; break; + } +} + +//------------------------------------------------------------------------ +/** + *

Directories which are searched for the musrview web assets, in order. + * + * \return the list of candidate directories + */ +std::vector WebDirCandidates() +{ + std::vector candidates; + + // 1. explicit override, e.g. when testing a modified page + if (const char *env = gSystem->Getenv("MUSRFIT_WEB_DIR")) + candidates.emplace_back(env); + + // 2. where 'cmake --install' put the page +#ifdef MUSRVIEW_WEB_DIR + candidates.emplace_back(MUSRVIEW_WEB_DIR); +#endif + + // 3. the source tree, so that a freshly built musrview also runs straight + // from the build directory, i.e. before anything has been installed +#ifdef MUSRVIEW_WEB_DIR_DEV + candidates.emplace_back(MUSRVIEW_WEB_DIR_DEV); +#endif + + // 4. common installation prefixes + if (const char *root = gSystem->Getenv("ROOTSYS")) + candidates.emplace_back(std::string(root) + "/share/musrfit/musrview"); + + candidates.emplace_back("/usr/local/share/musrfit/musrview"); + candidates.emplace_back("/usr/share/musrfit/musrview"); + + return candidates; +} + +//------------------------------------------------------------------------ +/** + *

Locates the directory holding the musrview web assets (the html page). + * + * \return the directory, or an empty string if none was found + */ +std::string FindWebDir() +{ + for (const auto &dir : WebDirCandidates()) { + if (dir.empty()) + continue; + const std::string page = dir + "/musrview.html"; + if (gSystem->AccessPathName(page.c_str()) == kFALSE) + return dir; + } + + return std::string(); +} + +} // anonymous namespace + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Creates the web front end of one PLOT block. + * + * \param plotNumber number of the PLOT block + * \param title msr-file title + * \param ww window width in pixels + * \param wh window height in pixels + * \param fourierDefault Fourier settings of the startup file + * \param markerList list of markers of the startup file + * \param colorList list of colors of the startup file + * \param batch true if no window shall be opened + * \param fourier true to start with the Fourier representation + * \param avg true to start with the averaged representation + * \param theoAsData true to calculate the theory only at the data points + */ +PMusrWebCanvas::PMusrWebCanvas(Int_t plotNumber, const std::string &title, Int_t ww, Int_t wh, + PMsrFourierStructure fourierDefault, + const PIntVector &markerList, const PIntVector &colorList, + Bool_t batch, Bool_t fourier, Bool_t avg, Bool_t theoAsData) : + fPlotNumber(plotNumber), fTitle(title), fWidth(ww), fHeight(wh), fBatchMode(batch), + fStartWithFourier(fourier), fStartWithAvg(avg) +{ + fModel = std::make_unique(plotNumber, fourierDefault, markerList, colorList, theoAsData); + + // musrview never wants the ROOT statistics/title boxes + gStyle->SetOptStat(0); + gStyle->SetOptTitle(0); + gStyle->SetHistMinimumZero(kTRUE); // needed to enforce proper bar option handling +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Releases the window and the drawing objects. The histograms themselves are + * owned by the model. + */ +PMusrWebCanvas::~PMusrWebCanvas() +{ + Close(); +} + +//-------------------------------------------------------------------------- +// Close (public) +//-------------------------------------------------------------------------- +/** + *

Closes the window and releases it. + * + *

This is not optional: as long as an RWebWindow exists, its http server and + * - for the embedded display - the Qt event loop keep running, and musrview + * would never return from main(). + */ +void PMusrWebCanvas::Close() +{ + if (!fWindow) + return; + + // do not let the disconnect callback fire into a half destroyed object + fWindow->SetDisconnectCallBack(nullptr); + fWindow->SetDataCallBack(nullptr); + fWindow->CloseConnections(); + fWindow.reset(); + + fDone = true; +} + +//-------------------------------------------------------------------------- +// SetMsrHandler (public) +//-------------------------------------------------------------------------- +/** + * \param msrHandler pointer to the msr-file handler + */ +void PMusrWebCanvas::SetMsrHandler(PMsrHandler *msrHandler) +{ + fModel->SetMsrHandler(msrHandler); +} + +//-------------------------------------------------------------------------- +// SetRunListCollection (public) +//-------------------------------------------------------------------------- +/** + * \param runList pointer to the run collection + */ +void PMusrWebCanvas::SetRunListCollection(PRunListCollection *runList) +{ + fModel->SetRunListCollection(runList); +} + +//-------------------------------------------------------------------------- +// Build (public) +//-------------------------------------------------------------------------- +/** + *

Builds the data of this plot block and derives the initial view mode. + * + * \return true if everything could be built + */ +Bool_t PMusrWebCanvas::Build() +{ + if (!fModel->BuildData()) { + fLastError = fModel->GetLastError(); + fValid = false; + return false; + } + + fMode = fModel->CreateDefaultMode(fStartWithFourier, fStartWithAvg); + + if (!fModel->Prepare(fMode)) { + fLastError = fModel->GetLastError(); + fValid = false; + return false; + } + + if (fTitle.empty()) + fTitle = fModel->GetTitle(); + + fValid = true; + + return true; +} + +//-------------------------------------------------------------------------- +// CollectDrawables (private) +//-------------------------------------------------------------------------- +/** + *

Builds the ordered list of objects which the current view mode shows, + * each with the ROOT draw option it has to be drawn with. + * + *

Which histograms these are is entirely decided by PMusrViewModel, this + * routine only knows how to present a set of curves. That is the reason why it + * is this short compared to the six Plot*() routines it replaces. Both the web + * front end and the graphics export go through it, so the two can never show + * different things. + * + * \param out where the (object, draw option) pairs are appended + * + * \return true on success + */ +Bool_t PMusrWebCanvas::CollectDrawables(std::vector> &out) +{ + out.clear(); + + if (!fValid) + return false; + + fFrame.reset(); + fMultiGraph.reset(); + fPlotLegend.reset(); + fPhaseLatex.reset(); + fRRFLatex.reset(); + + std::string xTitle, yTitle; + fModel->GetAxisTitles(fMode, xTitle, yTitle); + + Double_t xmin = 0.0, xmax = 0.0, ymin = 0.0, ymax = 0.0; + fModel->GetDefaultXRange(fMode, xmin, xmax); + fModel->GetDefaultYRange(fMode, xmin, xmax, ymin, ymax); + + // logarithmic axes are pad properties; for the web front end they travel as + // draw options of the first object, which is what JSROOT expects + TString padOpt(""); + if (fModel->IsLogX()) + padOpt += "logx"; + if (fModel->IsLogY()) + padOpt += "logy"; + + if (fModel->IsNonMusr()) { + //-- non-muSR: a TMultiGraph of error graphs --------------------------- + std::vector data = fModel->GetNonMusrCurves(fMode, kFALSE); + std::vector theo = fModel->GetNonMusrCurves(fMode, kTRUE); + + if (data.empty()) + return false; + + fMultiGraph = std::make_unique(); + for (auto *g : data) { + // TMultiGraph takes ownership, hence hand it a copy + TGraphErrors *ge = new TGraphErrors(*g); + ge->SetEditable(false); + fMultiGraph->Add(ge, "p"); + } + for (auto *g : theo) { + TGraphErrors *ge = new TGraphErrors(*g); + ge->SetEditable(false); + fMultiGraph->Add(ge, "l"); + } + // A TMultiGraph only gets its axes when it is drawn, so the ranges and the + // axis titles have to be put onto its frame histogram here. GetHistogram() + // creates that frame without needing a pad. + TH1F *mgFrame = fMultiGraph->GetHistogram(); + if (mgFrame) { + mgFrame->SetMinimum(ymin); + mgFrame->SetMaximum(ymax); + mgFrame->GetXaxis()->SetLimits(xmin, xmax); + // with more than one data set the axes cannot be labelled unambiguously, + // there a legend names what is plotted against what + if (data.size() == 1) { + mgFrame->GetXaxis()->SetTitle(xTitle.c_str()); + mgFrame->GetYaxis()->SetTitle(yTitle.c_str()); + mgFrame->GetYaxis()->SetTitleOffset(1.3); + } + } + + out.emplace_back(fMultiGraph.get(), TString("a") + padOpt); + + const std::vector &labels = fModel->GetInfoLines(); + if (data.size() > 1) { + fPlotLegend = std::make_unique(0.8, 0.8, 1.0, 1.0); + for (size_t i=0; iAddEntry(data[i], i < labels.size() ? labels[i].c_str() : "", "p"); + out.emplace_back(fPlotLegend.get(), TString("")); + } + + return true; + } + + //-- muSR: histograms on a common frame --------------------------------- + const PMusrViewCurves curves = fModel->GetCurves(fMode); + if (!curves.fValid) { + fLastError = "the requested view is not available."; + return false; + } + + // the frame carries the axes; the individual histograms may have different + // binning and ranges, hence they are drawn on top of it with 'same' + UInt_t noOfPoints = 1000; + for (auto *h : curves.fData) { + if (h && (h->GetNbinsX() > static_cast(noOfPoints))) + noOfPoints = h->GetNbinsX(); + } + noOfPoints *= 2; + + fFrame = std::make_unique("musrviewFrame", "", noOfPoints, xmin, xmax); + fFrame->SetDirectory(nullptr); + fFrame->SetStats(kFALSE); + fFrame->SetMinimum(ymin); + fFrame->SetMaximum(ymax); + fFrame->GetYaxis()->SetLimits(ymin, ymax); + fFrame->GetXaxis()->SetTitle(xTitle.c_str()); + fFrame->GetYaxis()->SetTitle(yTitle.c_str()); + fFrame->GetYaxis()->SetTitleOffset(1.3); + // an empty option, as TPad::DrawFrame() uses, so that the graphics export + // keeps producing exactly the same picture as before + out.emplace_back(fFrame.get(), padOpt); + + // Drawing options per representation, kept as they have always been: + // - time domain spectra carry meaningful per-bin errors, hence "pe" + // - Fourier spectra and averaged spectra do not: a freshly filled TH1F + // reports sqrt(content) as bin error, which would drown the plot in + // meaningless error bars, hence plain "p" + // - the theory is a line in the time domain and a histogram outline in the + // frequency domain + TString dataOpt("pesame"); + TString theoOpt("lsame"); + if (fMode.IsAverage()) { + dataOpt = "psame"; + theoOpt = "same"; + } else if (fMode.IsFourier()) { + dataOpt = fMode.IsDifference() ? "plsame" : "psame"; + theoOpt = "same"; + } + + for (size_t i=0; iGetRRFText(); + if (!rrf.empty()) { + fRRFLatex = std::make_unique(0.1, 0.92, rrf.c_str()); + fRRFLatex->SetNDC(kTRUE); + fRRFLatex->SetTextFont(62); + fRRFLatex->SetTextSize(0.03); + out.emplace_back(fRRFLatex.get(), TString("same")); + } + + // Fourier phase label; it describes the phase of the individual spectra and + // is therefore not shown in the averaged view + const std::string phase = fMode.IsAverage() ? std::string() : fModel->GetFourierPhaseText(fMode); + if (!phase.empty()) { + fPhaseLatex = std::make_unique(0.7, 0.85, phase.c_str()); + fPhaseLatex->SetNDC(kTRUE); + fPhaseLatex->SetTextFont(62); + fPhaseLatex->SetTextSize(0.03); + out.emplace_back(fPhaseLatex.get(), TString("same")); + } + + return true; +} + +//-------------------------------------------------------------------------- +// FillPlotPad (private) +//-------------------------------------------------------------------------- +/** + *

Draws the spectra of the current view mode into the given pad. Used by the + * graphics export, which runs with ROOT in batch mode. + * + * \param pad the pad to draw into + * + * \return true on success + */ +Bool_t PMusrWebCanvas::FillPlotPad(TVirtualPad *pad) +{ + if (pad == nullptr) + return false; + + std::vector> drawables; + if (!CollectDrawables(drawables)) + return false; + + pad->cd(); + pad->SetLogx(fModel->IsLogX() ? 1 : 0); + pad->SetLogy(fModel->IsLogY() ? 1 : 0); + + for (auto &entry : drawables) { + // the log flags are pad properties here, not draw options + TString opt(entry.second); + opt.ReplaceAll("logx", ""); + opt.ReplaceAll("logy", ""); + entry.first->Draw(opt.Data()); + } + + return true; +} + +//-------------------------------------------------------------------------- +// BuildPlotJson (private) +//-------------------------------------------------------------------------- +/** + *

Serializes the objects of the current view as a TList, together with the + * draw option of each of them. + * + *

Deliberately no TCanvas is built here. A TCanvas creates and maps a window + * in its constructor, which would put an empty ROOT canvas next to the musrview + * window, and the ways around that either crash (toggling gROOT batch state) or + * open yet another window (ROOT web canvases). Since JSROOT does the painting + * anyway, sending the drawables plus their options is both simpler and closer + * to the ROOT7 idea of decoupling data from its graphics representation. + * + * \param json where the JSON is returned + * + * \return true on success + */ +Bool_t PMusrWebCanvas::BuildPlotJson(TString &json) +{ + std::vector> drawables; + if (!CollectDrawables(drawables)) + return false; + + // a non-owning list: everything in it belongs to the model or to this object + TList list; + for (auto &entry : drawables) + list.Add(entry.first, entry.second.Data()); + + json = TBufferJSON::ToJSON(&list, TBufferJSON::kNoSpaces + TBufferJSON::kSameSuppression); + + list.Clear("nodelete"); + + return !json.IsNull(); +} + +//-------------------------------------------------------------------------- +// JsonEscape (static private) +//-------------------------------------------------------------------------- +/** + *

Escapes a string so that it can be used as a JSON string value. + * + * \param str the string to escape + * + * \return the escaped string + */ +std::string PMusrWebCanvas::JsonEscape(const std::string &str) +{ + std::string out; + out.reserve(str.size()+8); + + for (unsigned char c : str) { + switch (c) { + case '"': out += "\\\""; break; + case '\\': out += "\\\\"; break; + case '\b': out += "\\b"; break; + case '\f': out += "\\f"; break; + case '\n': out += "\\n"; break; + case '\r': out += "\\r"; break; + case '\t': out += "\\t"; break; + default: + if (c < 0x20) { + char buf[8]; + snprintf(buf, sizeof(buf), "\\u%04x", c); + out += buf; + } else { + out += static_cast(c); + } + break; + } + } + + return out; +} + +//-------------------------------------------------------------------------- +// JsonArray (static private) +//-------------------------------------------------------------------------- +/** + *

Turns a vector of strings into a JSON array. + * + * \param vec the strings + * + * \return the JSON array as string + */ +std::string PMusrWebCanvas::JsonArray(const std::vector &vec) +{ + std::string out = "["; + for (size_t i=0; iAssembles everything the client needs besides the spectra: the view mode + * (so that the toolbar can render itself without duplicating any logic) and the + * text blocks of the panels. + * + * \return the JSON object as string + */ +std::string PMusrWebCanvas::BuildStateJson() const +{ + std::string json = "{"; + json += "\"plotNumber\":" + std::to_string(fPlotNumber) + ","; + json += "\"title\":\"" + JsonEscape(fTitle) + "\","; + json += "\"mode\":" + fMode.ToJson() + ","; + json += "\"statistic\":\"" + JsonEscape(fModel->GetStatisticLine()) + "\","; + json += "\"rrf\":\"" + JsonEscape(fModel->GetRRFText()) + "\","; + json += "\"parameters\":" + JsonArray(fModel->GetParameterLines()) + ","; + json += "\"theory\":" + JsonArray(fModel->GetTheoryLines()) + ","; + json += "\"info\":" + JsonArray(fModel->GetInfoLines()) + ","; + json += "\"runs\":" + BuildRunsJson() + ","; + json += "\"error\":\"" + JsonEscape(fLastError) + "\""; + json += "}"; + + return json; +} + +//-------------------------------------------------------------------------- +// BuildRunsJson (private) +//-------------------------------------------------------------------------- +/** + *

Describes every run of the plot block by its info line and by the symbol + * it is plotted with. + * + *

The former TLegend drew those symbols itself. The run info is an HTML + * panel now, so colour and marker shape have to travel to the page, which draws + * a matching symbol next to each line. + * + * \return the JSON array as string + */ +std::string PMusrWebCanvas::BuildRunsJson() const +{ + const std::vector &labels = fModel->GetInfoLines(); + + std::vector> style; // (marker style, marker colour) + if (fModel->IsNonMusr()) { + for (const auto &ds : fModel->GetNonMusrData()) + if (ds.data) + style.emplace_back(ds.data->GetMarkerStyle(), ds.data->GetMarkerColor()); + } else { + for (const auto &ds : fModel->GetData()) + if (ds.data) + style.emplace_back(ds.data->GetMarkerStyle(), ds.data->GetMarkerColor()); + } + + std::string json = "["; + for (size_t i=0; iSends the view state and the panel contents to the client. + * + * \param connid connection the message is sent to + */ +void PMusrWebCanvas::SendState(unsigned connid) +{ + if (!fWindow) + return; + const std::string payload = "STATE:" + BuildStateJson(); + if (gSystem->Getenv("MUSRVIEW_DEBUG")) + std::cerr << ">> musrview DEBUG send STATE (" << payload.size() << " bytes)" << std::endl; + fWindow->Send(connid, payload); +} + +//-------------------------------------------------------------------------- +// SendPlot (private) +//-------------------------------------------------------------------------- +/** + *

Serializes the spectra canvas with TBufferJSON and sends it to the client, + * where JSROOT renders it. + * + * \param connid connection the message is sent to + */ +void PMusrWebCanvas::SendPlot(unsigned connid) +{ + if (!fWindow) + return; + + // Compact level 23 = kNoSpaces + kSameSuppression. The type names must be + // kept (i.e. no kSkipTypeInfo), otherwise JSROOT cannot tell what it got. + TString json; + if (!BuildPlotJson(json)) { + fWindow->Send(connid, "ERROR:" + JsonEscape(fLastError)); + return; + } + + if (gSystem->Getenv("MUSRVIEW_DEBUG")) + std::cerr << ">> musrview DEBUG send PLOT (" << json.Length() << " bytes), " + << gROOT->GetListOfCanvases()->GetSize() << " ROOT canvas(es) alive" << std::endl; + fWindow->Send(connid, "PLOT:" + std::string(json.Data())); +} + +//-------------------------------------------------------------------------- +// SendAll (private) +//-------------------------------------------------------------------------- +/** + *

Sends state and spectra in one go. + * + * \param connid connection the messages are sent to + */ +void PMusrWebCanvas::SendAll(unsigned connid) +{ + SendState(connid); + SendPlot(connid); +} + +//-------------------------------------------------------------------------- +// ApplyCommand (private) +//-------------------------------------------------------------------------- +/** + *

Applies a command coming from the page onto the view mode. + * + *

Every command is a single, well defined operation on PMusrViewMode. The + * combinatorial explosion which used to live in HandleCmdKey()/HandleMenuPopup() + * is gone: the mode object itself knows which combinations exist and the model + * knows how to produce them. + * + * \param cmd the command + * + * \return true if the spectra have to be redrawn + */ +Bool_t PMusrWebCanvas::ApplyCommand(const std::string &cmd) +{ + Bool_t changed = kFALSE; + + if (cmd == "difference") { + changed = fMode.ToggleDifference(); + } else if (cmd == "fourier") { + changed = fMode.ToggleFourier(); + } else if (cmd == "average") { + changed = fMode.ToggleAverage(); + if (!changed && !fMode.IsAverageAllowed()) + fLastError = "averaging a single data set does not make sense."; + } else if (cmd == "data") { + changed = fMode.SetDomain(EDomain::kTime); + } else if (cmd == "fourier-real") { + changed = fMode.SetFourierComponent(EFourierComponent::kReal); + } else if (cmd == "fourier-imag") { + changed = fMode.SetFourierComponent(EFourierComponent::kImag); + } else if (cmd == "fourier-real-imag") { + changed = fMode.SetFourierComponent(EFourierComponent::kRealAndImag); + } else if (cmd == "fourier-power") { + changed = fMode.SetFourierComponent(EFourierComponent::kPower); + } else if (cmd == "fourier-phase") { + changed = fMode.SetFourierComponent(EFourierComponent::kPhase); + } else if (cmd == "fourier-phase-opt-real") { + changed = fMode.SetFourierComponent(EFourierComponent::kPhaseOptReal); + } else if (cmd == "phase+") { + fModel->IncrementFourierPhase(fMode); + changed = fMode.IsPhaseAdjustable(); + } else if (cmd == "phase-") { + fModel->DecrementFourierPhase(fMode); + changed = fMode.IsPhaseAdjustable(); + } else if (cmd == "unzoom") { + changed = kTRUE; + } else { + fLastError = "unknown command '" + cmd + "'."; + return kFALSE; + } + + if (changed) { + fLastError.clear(); + if (!fModel->Prepare(fMode)) + fLastError = fModel->GetLastError(); + } + + return changed; +} + +//-------------------------------------------------------------------------- +// ProcessMessage (private) +//-------------------------------------------------------------------------- +/** + *

Handles a message coming from the page. + * + * \param connid connection the message came from + * \param msg the message + */ +void PMusrWebCanvas::ProcessMessage(unsigned connid, const std::string &msg) +{ + if (gSystem->Getenv("MUSRVIEW_DEBUG")) + std::cerr << ">> musrview DEBUG recv[" << connid << "]: " << msg.substr(0,80) << std::endl; + + if (msg == "init") { + SendAll(connid); + return; + } + + if (msg == "quit") { + fDone = true; + if (fDoneCallback) + fDoneCallback(0); + return; + } + + if (msg.rfind("cmd:", 0) == 0) { + const std::string cmd = msg.substr(4); + const Bool_t changed = ApplyCommand(cmd); + SendState(connid); + if (changed) + SendPlot(connid); + return; + } + + if (msg.rfind("zoom:", 0) == 0) { + // zoom:: - the Fourier transform follows the visible time range + const std::string arg = msg.substr(5); + const size_t sep = arg.find(':'); + if (sep == std::string::npos) + return; + try { + const Double_t xmin = std::stod(arg.substr(0, sep)); + const Double_t xmax = std::stod(arg.substr(sep+1)); + if (fMode.IsTimeDomain()) + fModel->SetVisibleTimeRange(xmin, xmax); + } catch (const std::exception &) { + // malformed message, ignore + } + return; + } + + if (msg.rfind("export:", 0) == 0) { + const std::string fileName = msg.substr(7); + if (ExportData(fileName)) + fWindow->Send(connid, "INFO:data exported to " + JsonEscape(fileName)); + else + fWindow->Send(connid, "ERROR:" + JsonEscape(fLastError)); + return; + } + + if (msg.rfind("save:", 0) == 0) { + const std::string arg = msg.substr(5); + const size_t sep = arg.find(':'); + if (sep == std::string::npos) + return; + const std::string fileName = arg.substr(0, sep); + const std::string format = arg.substr(sep+1); + if (SaveGraphics(fileName, format)) + fWindow->Send(connid, "INFO:graphics written (" + JsonEscape(format) + ")"); + else + fWindow->Send(connid, "ERROR:" + JsonEscape(fLastError)); + return; + } +} + +//-------------------------------------------------------------------------- +// Show (public) +//-------------------------------------------------------------------------- +/** + *

Opens the window. + * + * \param where display to be used: an empty string picks the user's standard + * browser, "qt6" opens an embedded Qt6 WebEngine window so that musrview + * appears as a stand-alone application, further values such as "chrome", + * "firefox" or "server" are passed on to ROOT unchanged. + * + * \return true if the window could be opened + */ +Bool_t PMusrWebCanvas::Show(const std::string &where) +{ + if (!fValid) { + fLastError = "canvas is not valid, refusing to show it."; + return false; + } + + const std::string webDir = FindWebDir(); + if (webDir.empty()) { + std::ostringstream err; + err << "couldn't locate the musrview web assets (musrview.html). Looked in:"; + for (const auto &dir : WebDirCandidates()) + if (!dir.empty()) + err << "\n>> " << dir; + err << "\n>> Set MUSRFIT_WEB_DIR to the directory holding musrview.html."; + fLastError = err.str(); + std::cerr << std::endl << ">> PMusrWebCanvas::Show(): **ERROR** " << fLastError << std::endl; + return false; + } + + fWindow = ROOT::RWebWindow::Create(); + if (!fWindow) { + fLastError = "couldn't create the web window."; + return false; + } + + fWindow->SetDefaultPage("file:" + webDir + "/musrview.html"); + fWindow->SetGeometry(fWidth, fHeight); + fWindow->SetConnLimit(1); + + fWindow->SetDataCallBack([this](unsigned connid, const std::string &msg) { + ProcessMessage(connid, msg); + }); + + fWindow->SetConnectCallBack([this](unsigned connid) { + if (gSystem->Getenv("MUSRVIEW_DEBUG")) + std::cerr << ">> musrview DEBUG connected: " << connid << std::endl; + }); + + fWindow->SetDisconnectCallBack([this](unsigned) { + fDone = true; + if (fDoneCallback) + fDoneCallback(0); + }); + + fWindow->Show(where); + + return true; +} + +//-------------------------------------------------------------------------- +// SaveGraphics (public) +//-------------------------------------------------------------------------- +/** + *

Writes the classic musrview layout (title, spectra, fit parameters, + * theory, run info) into a graphics file. + * + *

This path deliberately does not involve the browser: batch operation must + * not depend on a browser installation, and the output stays reproducible. + * + * \param msrFileName name of the msr-file; the output name is derived from it + * \param graphicsFormat extension such as "png", "pdf", "eps", ... + * + * \return true on success + */ +Bool_t PMusrWebCanvas::SaveGraphics(const std::string &msrFileName, const std::string &graphicsFormat) +{ + if (!fValid) { + fLastError = "canvas is not valid."; + return false; + } + + // derive the output file name: _[_F]. + TString str(msrFileName.c_str()); + Int_t idx = -1; + Int_t size = 0; + if (str.Contains(".msr")) { + idx = str.Index(".msr"); + size = 4; + } else if (str.Contains(".mlog")) { + idx = str.Index(".mlog"); + size = 5; + } + if (idx == -1) { + fLastError = "file name '" + msrFileName + "' is invalid."; + std::cerr << std::endl << ">> PMusrWebCanvas::SaveGraphics(): **ERROR** " << fLastError << std::endl; + return false; + } + + Char_t ext[64]; + if (fMode.IsFourier()) + snprintf(ext, sizeof(ext), "_%d_F", fPlotNumber); + else + snprintf(ext, sizeof(ext), "_%d", fPlotNumber); + str.Replace(idx, size, ext, strlen(ext)); + idx += strlen(ext); + size = strlen(ext); + snprintf(ext, sizeof(ext), ".%s", graphicsFormat.c_str()); + str.Replace(idx, size, ext, strlen(ext)); + + // Build the classic layout. Its size is deliberately not the size of the + // musrview window: the exported graphics have always been 800x600 and users + // (and the reference images of the test suite) rely on that. + const Int_t exportWidth = 800; + const Int_t exportHeight = 600; + + std::string name = "musrviewBatch" + std::to_string(fPlotNumber); + TCanvas canvas(name.c_str(), fTitle.c_str(), exportWidth, exportHeight); + canvas.SetBatch(kTRUE); + + const Double_t yInfo = 0.1; + const Double_t yTitle = 0.95; + const Double_t xTheo = 0.75; + const Int_t white = TColor::GetColor(255, 255, 255); + + TPaveText titlePad(0.0, yTitle, 1.0, 1.0, "NDC"); + titlePad.SetFillColor(white); + titlePad.SetTextAlign(12); // middle, left + titlePad.AddText(fTitle.c_str()); + titlePad.Draw(); + + TPad dataTheoryPad("dataTheoryPad", "dataTheoryPad", 0.0, yInfo, xTheo, yTitle); + dataTheoryPad.SetFillColor(white); + dataTheoryPad.Draw(); + + TPaveText parameterPad(xTheo, 0.5, 1.0, yTitle, "NDC"); + parameterPad.SetFillColor(white); + parameterPad.SetTextAlign(13); // top, left + parameterPad.SetTextFont(102); // courier bold + { + const std::vector &lines = fModel->GetParameterLines(); + const Double_t yoffset = (lines.size() > 20) ? 1.0/(lines.size()+1) : 0.05; + for (size_t i=0; i &lines = fModel->GetTheoryLines(); + const Double_t yoffset = (lines.size() > 20) ? 1.0/(lines.size()+1) : 0.05; + for (size_t i=0; iGetStatisticLine().c_str()); + { + // the legend always lists the runs of the plot block, also in the averaged + // view where a single curve is shown + const std::vector &lines = fModel->GetInfoLines(); + const PMusrViewDataList &data = fModel->GetData(); + for (size_t i=0; i> musrview: writing " << str.Data() << std::endl; + canvas.SaveAs(str.Data()); + + // the primitives belong to the model/this function, do not let ROOT delete them + dataTheoryPad.GetListOfPrimitives()->Clear(); + canvas.GetListOfPrimitives()->Clear(); + + return true; +} + +//-------------------------------------------------------------------------- +// ExportData (public) +//-------------------------------------------------------------------------- +/** + *

Writes the currently shown data in ascii column format. + * + * \param fileName name of the file to be written + * + * \return true on success + */ +Bool_t PMusrWebCanvas::ExportData(const std::string &fileName) +{ + if (!fValid) { + fLastError = "canvas is not valid."; + return false; + } + + Double_t xmin = 0.0, xmax = 0.0; + fModel->GetDefaultXRange(fMode, xmin, xmax); + + // in the time domain the user may have zoomed in; export what is on screen + TH1F *frame = fModel->GetHistoFrame(); + if (fMode.IsTimeDomain() && frame) { + const Int_t first = frame->GetXaxis()->GetFirst(); + const Int_t last = frame->GetXaxis()->GetLast(); + xmin = frame->GetXaxis()->GetBinLowEdge(first); + xmax = frame->GetXaxis()->GetBinLowEdge(last) + frame->GetXaxis()->GetBinWidth(last); + } + + if (!fModel->ExportData(fileName.c_str(), fMode, xmin, xmax)) { + fLastError = fModel->GetLastError(); + return false; + } + + return true; +} + +} // namespace musrview diff --git a/src/include/PMusrViewMode.h b/src/include/PMusrViewMode.h new file mode 100644 index 00000000..ffb6e195 --- /dev/null +++ b/src/include/PMusrViewMode.h @@ -0,0 +1,199 @@ +/*************************************************************************** + + PMusrViewMode.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2026 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 _PMUSRVIEWMODE_H_ +#define _PMUSRVIEWMODE_H_ + +#include + +namespace musrview { + +//-------------------------------------------------------------------------- +/** + * \brief Representation domain of the currently shown spectra. + */ +enum class EDomain { + kTime, ///< time domain, i.e. the raw asymmetry/single histogram spectra + kFourier ///< frequency domain, i.e. the Fourier transform of the spectra +}; + +//-------------------------------------------------------------------------- +/** + * \brief Component of the complex Fourier transform which shall be shown. + * + * Only meaningful if the domain is EDomain::kFourier. + */ +enum class EFourierComponent { + kReal, ///< real part + kImag, ///< imaginary part + kRealAndImag, ///< real and imaginary part in a single plot + kPower, ///< power spectrum, i.e. |F|^2 + kPhase, ///< phase spectrum, i.e. arg(F) + kPhaseOptReal ///< real part with an optimized phase +}; + +//-------------------------------------------------------------------------- +/** + * \brief Position of the difference operation relative to the Fourier transform. + * + * The difference is always 'data - theory'. In the time domain the placement is + * irrelevant, hence kBeforeFourier is used as the canonical value there. + */ +enum class EDiffStage { + kNone, ///< no difference is formed at all + kBeforeFourier, ///< difference is formed first, the result is Fourier transformed (d-f) + kAfterFourier ///< data and theory are Fourier transformed first, then subtracted (f-d) +}; + +//-------------------------------------------------------------------------- +/** + * \brief Fully decoupled description of what musrview currently shows. + * + *

The plot view of musrview is spanned by four *independent* properties: + * + * + * + * + * + * + * + *
propertyvalues
domaintime / Fourier
Fourier componentreal / imag / real+imag / power / phase / phase-opt-real
differenceoff / on, and if on: formed before or after the Fourier transform
averageoff / on
+ * + *

Previously these properties were encoded in a set of entangled variables + * (\c fCurrentPlotView, \c fDifferenceView, \c fAveragedView plus a \em static + * \c lastKeySwitch which silently leaked state between canvases). Every new + * combination required another branch in an already deeply nested cascade, + * which is why combinations like 'Fourier-Difference-Average' were fragile. + * + *

This class keeps the four properties orthogonal. The only coupling which + * genuinely exists in the physics is expressed explicitly: + * + * - the Fourier component only matters in the Fourier domain, + * - the difference stage only matters if difference *and* Fourier are active, + * - a phase can only be adjusted for the real/imaginary components, + * - some plot types (non-muSR) cannot be Fourier transformed at all, and + * averaging requires more than one data set. + * + *

The latter two are capabilities of the underlying data, not of the view, + * and are therefore injected via SetFourierAllowed()/SetAverageAllowed(). All + * toggles are no-ops when the corresponding capability is missing, so the + * caller never has to guard them. + * + *

The class is deliberately free of any ROOT and GUI dependency so that it + * can be used by the web front end, the batch/graphics export and the tests + * alike. + */ +class PMusrViewMode +{ + public: + PMusrViewMode() = default; + + //-- capabilities of the underlying data ------------------------------- + /// Allow/forbid the Fourier domain (forbidden for non-muSR plots). + void SetFourierAllowed(bool allowed); + /// Allow/forbid averaging (needs more than one data set). + void SetAverageAllowed(bool allowed); + /// Returns true if the Fourier domain can be entered. + bool IsFourierAllowed() const { return fFourierAllowed; } + /// Returns true if averaging is possible. + bool IsAverageAllowed() const { return fAverageAllowed; } + + //-- state modifiers --------------------------------------------------- + /// Switch between time and Fourier domain ('f'). Returns true if the state changed. + bool ToggleFourier(); + /// Switch the difference on/off ('d'). Returns true if the state changed. + bool ToggleDifference(); + /// Switch averaging on/off ('a'). Returns true if the state changed. + bool ToggleAverage(); + + /// Explicitly select the domain. Returns true if the state changed. + bool SetDomain(EDomain domain); + /// Explicitly switch the difference on/off. Returns true if the state changed. + bool SetDifference(bool on); + /// Explicitly switch averaging on/off. Returns true if the state changed. + bool SetAverage(bool on); + /// Select the Fourier component; implies the Fourier domain. Returns true if the state changed. + bool SetFourierComponent(EFourierComponent comp); + + /// Reset to the plain time domain view, keeping the capabilities. + void Reset(); + + //-- state queries ----------------------------------------------------- + EDomain GetDomain() const { return fDomain; } + EFourierComponent GetFourierComponent() const { return fFourierComponent; } + EDiffStage GetDiffStage() const; + + bool IsTimeDomain() const { return fDomain == EDomain::kTime; } + bool IsFourier() const { return fDomain == EDomain::kFourier; } + bool IsDifference() const { return fDifference; } + bool IsAverage() const { return fAverage; } + + /// True if the difference is formed in the time domain and then transformed (d-f). + bool IsDiffThenFourier() const { return GetDiffStage() == EDiffStage::kBeforeFourier && IsFourier(); } + /// True if data and theory are transformed and then subtracted (f-d). + bool IsFourierThenDiff() const { return GetDiffStage() == EDiffStage::kAfterFourier; } + + /// True if a Fourier phase can be applied in the current state. + bool IsPhaseAdjustable() const; + + //-- presentation ------------------------------------------------------ + /// Stable, machine readable key, e.g. "fourier-real+diff-fd+avg". + std::string GetKey() const; + /// Human readable description, e.g. "Fourier real part, difference (f-d), averaged". + std::string GetLabel() const; + /// JSON object describing the state, consumed by the web front end. + std::string ToJson() const; + + /// Name of a Fourier component, e.g. "real". + static std::string GetFourierComponentName(EFourierComponent comp); + /// Human readable name of a Fourier component, e.g. "real part". + static std::string GetFourierComponentLabel(EFourierComponent comp); + /// Maps the FOURIER_PLOT_* tag of the msr-file onto a Fourier component. + static EFourierComponent FourierComponentFromMsrTag(int plotTag); + /// Maps a Fourier component onto the FOURIER_PLOT_* tag of the msr-file. + static int MsrTagFromFourierComponent(EFourierComponent comp); + + /// Equality, capabilities are deliberately not part of the comparison. + bool operator==(const PMusrViewMode &rhs) const; + bool operator!=(const PMusrViewMode &rhs) const { return !(*this == rhs); } + + private: + EDomain fDomain{EDomain::kTime}; ///< time or Fourier domain + EFourierComponent fFourierComponent{EFourierComponent::kPower}; ///< which part of the FT is shown + bool fDifference{false}; ///< is the difference data-theory shown? + bool fDiffBeforeFourier{true}; ///< d-f (true) or f-d (false); only used if difference && Fourier + bool fAverage{false}; ///< is the average over all data sets shown? + + bool fFourierAllowed{true}; ///< can the Fourier domain be entered at all? + bool fAverageAllowed{true}; ///< can the data be averaged at all? +}; + +} // namespace musrview + +#endif // _PMUSRVIEWMODE_H_ diff --git a/src/include/PMusrViewModel.h b/src/include/PMusrViewModel.h new file mode 100644 index 00000000..3111042e --- /dev/null +++ b/src/include/PMusrViewModel.h @@ -0,0 +1,391 @@ +/*************************************************************************** + + PMusrViewModel.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2026 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 _PMUSRVIEWMODEL_H_ +#define _PMUSRVIEWMODEL_H_ + +#include +#include + +#include +#include +#include + +#include "PMusr.h" +#include "PMusrViewMode.h" +#ifndef __MAKECLING__ +#include "PMsrHandler.h" +#include "PRunListCollection.h" +#endif // __MAKECLING__ + +namespace musrview { + +//------------------------------------------------------------------------ +/** + * \brief Keeps an explicitly given plot range of the msr-file PLOT block. + */ +class PMusrViewPlotRange +{ + public: + PMusrViewPlotRange() = default; + + void SetXRange(Double_t xmin, Double_t xmax); + void SetYRange(Double_t ymin, Double_t ymax); + + Bool_t IsXRangePresent() const { return fXRangePresent; } + Bool_t IsYRangePresent() const { return fYRangePresent; } + + Double_t GetXmin() const { return fXmin; } + Double_t GetXmax() const { return fXmax; } + Double_t GetYmin() const { return fYmin; } + Double_t GetYmax() const { return fYmax; } + + private: + Bool_t fXRangePresent{false}; ///< flag: true if x-range explicitly set + Bool_t fYRangePresent{false}; ///< flag: true if y-range explicitly set + Double_t fXmin{0.0}; ///< minimum x value + Double_t fXmax{0.0}; ///< maximum x value + Double_t fYmin{0.0}; ///< minimum y value + Double_t fYmax{0.0}; ///< maximum y value +}; + +//------------------------------------------------------------------------ +/** + *

All histograms belonging to a single entry of the PLOT block for the fit + * types single histogram and asymmetry. + */ +struct PMusrViewDataSet { + TH1F *data{nullptr}; ///< data histogram + TH1F *dataFourierRe{nullptr}; ///< real part of the Fourier transform of the data histogram + TH1F *dataFourierIm{nullptr}; ///< imaginary part of the Fourier transform of the data histogram + TH1F *dataFourierPwr{nullptr}; ///< power spectrum of the Fourier transform of the data histogram + TH1F *dataFourierPhase{nullptr}; ///< phase spectrum of the Fourier transform of the data histogram + TH1F *dataFourierPhaseOptReal{nullptr}; ///< phase optimized real part spectrum of the data histogram + TH1F *theory{nullptr}; ///< theory histogram belonging to the data histogram + TH1F *theoryFourierRe{nullptr}; ///< real part of the Fourier transform of the theory histogram + TH1F *theoryFourierIm{nullptr}; ///< imaginary part of the Fourier transform of the theory histogram + TH1F *theoryFourierPwr{nullptr}; ///< power spectrum of the Fourier transform of the theory histogram + TH1F *theoryFourierPhase{nullptr}; ///< phase spectrum of the Fourier transform of the theory histogram + TH1F *theoryFourierPhaseOptReal{nullptr}; ///< phase optimized real part spectrum of the theory histogram + TH1F *diff{nullptr}; ///< difference histogram, i.e. data-theory + TH1F *diffFourierRe{nullptr}; ///< real part of the Fourier transform of the diff histogram + TH1F *diffFourierIm{nullptr}; ///< imaginary part of the Fourier transform of the diff histogram + TH1F *diffFourierPwr{nullptr}; ///< power spectrum of the Fourier transform of the diff histogram + TH1F *diffFourierPhase{nullptr}; ///< phase spectrum of the Fourier transform of the diff histogram + TH1F *diffFourierPhaseOptReal{nullptr}; ///< phase optimized real part spectrum of the diff histogram + PMusrViewPlotRange *dataRange{nullptr}; ///< keep the msr-file plot data range + UInt_t diffFourierTag{0}; ///< 0=not relevant, 1=d-f (Fourier of difference), 2=f-d (difference of Fourier) +}; + +//------------------------------------------------------------------------ +/** + *

typedef to make the code more readable: list of histogram data sets. + */ +typedef std::vector PMusrViewDataList; + +//------------------------------------------------------------------------ +/** + *

All error graphs belonging to a single entry of the PLOT block for the fit + * type non-muSR. + */ +struct PMusrViewNonMusrDataSet { + TGraphErrors *data{nullptr}; ///< data error graph + TGraphErrors *dataFourierRe{nullptr}; ///< real part of the Fourier transform of the data error graph + TGraphErrors *dataFourierIm{nullptr}; ///< imaginary part of the Fourier transform of the data error graph + TGraphErrors *dataFourierPwr{nullptr}; ///< power spectrum of the Fourier transform of the data error graph + TGraphErrors *dataFourierPhase{nullptr}; ///< phase spectrum of the Fourier transform of the data error graph + TGraphErrors *theory{nullptr}; ///< theory graph belonging to the data error graph + TGraphErrors *theoryFourierRe{nullptr}; ///< real part of the Fourier transform of the theory error graph + TGraphErrors *theoryFourierIm{nullptr}; ///< imaginary part of the Fourier transform of the theory error graph + TGraphErrors *theoryFourierPwr{nullptr}; ///< power spectrum of the Fourier transform of the theory error graph + TGraphErrors *theoryFourierPhase{nullptr}; ///< phase spectrum of the Fourier transform of the theory error graph + TGraphErrors *diff{nullptr}; ///< difference error graph, i.e. data-theory + TGraphErrors *diffFourierRe{nullptr}; ///< real part of the Fourier transform of the diff error graph + TGraphErrors *diffFourierIm{nullptr}; ///< imaginary part of the Fourier transform of the diff error graph + TGraphErrors *diffFourierPwr{nullptr}; ///< power spectrum of the Fourier transform of the diff error graph + TGraphErrors *diffFourierPhase{nullptr}; ///< phase spectrum of the Fourier transform of the diff error graph + PMusrViewPlotRange *dataRange{nullptr}; ///< keep the msr-file plot data range + UInt_t diffFourierTag{0}; ///< 0=not relevant, 1=d-f, 2=f-d +}; + +//------------------------------------------------------------------------ +/** + *

typedef to make the code more readable: list of error graph data sets. + */ +typedef std::vector PMusrViewNonMusrDataList; + +//------------------------------------------------------------------------ +/** + *

Data structure needed for the ascii dump of musrview. + */ +struct PMusrViewAsciiDump { + PDoubleVector dataX; ///< x-axis data set + PDoubleVector data; ///< y-axis data set + PDoubleVector dataErr; ///< error of the y-axis data set +}; + +//------------------------------------------------------------------------ +/** + *

typedef to make the code more readable: vector of ascii dump data. + */ +typedef std::vector PMusrViewAsciiDumpVector; + +//------------------------------------------------------------------------ +/** + * \brief The set of curves which a given PMusrViewMode makes visible. + * + *

This is the single place which knows how a view mode maps onto the + * histograms of PMusrViewDataSet. Plotting, ascii export and graphics export + * all go through it, instead of each repeating the same mode cascade. + * + *

\c fDataSecond / \c fTheorySecond are only filled for the combined + * real+imaginary Fourier view, where two curves per data set are shown. + */ +struct PMusrViewCurves { + std::vector fData; ///< data curves to be shown + std::vector fTheory; ///< theory curves to be shown; empty in difference views + std::vector fDataSecond; ///< second data curve (imaginary part) for the real+imag view + std::vector fTheorySecond; ///< second theory curve (imaginary part) for the real+imag view + std::vector fLabel; ///< legend label per data curve + Bool_t fHasError{true}; ///< true if the y-errors of the data curves are meaningful + Bool_t fValid{false}; ///< true if the requested mode could be served + + /// Number of data curves. + size_t Size() const { return fData.size(); } + /// True if nothing is to be shown. + Bool_t IsEmpty() const { return fData.empty(); } +}; + +//-------------------------------------------------------------------------- +/** + * \brief GUI-free data layer of musrview. + * + *

PMusrViewModel owns everything which is shown by musrview and knows how to + * calculate it: the raw data/theory histograms, their Fourier transforms, the + * difference data-theory (in both orders, d-f and f-d), the average over all + * data sets, and the text blocks (title, fit parameters, theory, run info). + * + *

It deliberately contains no ROOT GUI code at all. The front end (web based + * in musrview, TCanvas based for the batch graphics export) asks the model for + * a PMusrViewCurves object matching the current PMusrViewMode and renders it. + * Adding a further representation therefore does not touch this class, and + * adding a further view mode does not touch any front end. + */ +class PMusrViewModel +{ + public: + PMusrViewModel(Int_t plotNumber, PMsrFourierStructure fourierDefault, + PIntVector markerList, PIntVector colorList, Bool_t theoAsData); + virtual ~PMusrViewModel(); + + PMusrViewModel(const PMusrViewModel &) = delete; + PMusrViewModel &operator=(const PMusrViewModel &) = delete; + +#ifndef __MAKECLING__ + /// Sets the msr-handler and takes over the FOURIER/RRF settings of the msr-file. + void SetMsrHandler(PMsrHandler *msrHandler); + /// Sets the run collection which delivers the fitted data. + void SetRunListCollection(PRunListCollection *runList) { fRunList = runList; } +#endif // __MAKECLING__ + + /// Builds all time domain histograms of the PLOT block. Returns false on error. + Bool_t BuildData(); + + /// True if the model is in a usable state. + Bool_t IsValid() const { return fValid; } + /// Last error message, empty if none. + const std::string &GetLastError() const { return fLastError; } + + /// MSR_PLOT_* tag of this plot block. + Int_t GetPlotType() const { return fPlotType; } + /// Plot block number. + Int_t GetPlotNumber() const { return fPlotNumber; } + /// True for the non-muSR fit type, which cannot be Fourier transformed. + Bool_t IsNonMusr() const { return fPlotType == MSR_PLOT_NON_MUSR; } + /// Number of data sets of this plot block. + size_t GetNoOfDataSets() const { return IsNonMusr() ? fNonMusrData.size() : fData.size(); } + + /// Creates a view mode which matches the capabilities of the loaded data. + PMusrViewMode CreateDefaultMode(Bool_t startWithFourier, Bool_t startWithAvg) const; + + /// Makes sure that everything the mode needs has been calculated. + Bool_t Prepare(const PMusrViewMode &mode); + + /// Tells the model which time range the user currently looks at. The + /// Fourier transform is calculated over exactly this window. + void SetVisibleTimeRange(Double_t xmin, Double_t xmax); + /// The frame histogram which defines the visible time range. + TH1F *GetHistoFrame() const { return fHistoFrame; } + + /// The curves which the given mode shows. Call Prepare() first. + PMusrViewCurves GetCurves(const PMusrViewMode &mode) const; + /// The non-muSR error graphs which the given mode shows. + std::vector GetNonMusrCurves(const PMusrViewMode &mode, Bool_t theory) const; + + /// Increments the Fourier phase by the configured increment. + void IncrementFourierPhase(const PMusrViewMode &mode); + /// Decrements the Fourier phase by the configured increment. + void DecrementFourierPhase(const PMusrViewMode &mode); + /// The currently applied Fourier phase(s). + const PDoubleVector &GetCurrentFourierPhase() const { return fCurrentFourierPhase; } + /// Text describing the current Fourier phase, empty if not applicable. + std::string GetFourierPhaseText(const PMusrViewMode &mode) const; + + /// x-axis range the given mode should show by default. + void GetDefaultXRange(const PMusrViewMode &mode, Double_t &xmin, Double_t &xmax) const; + /// y-axis range the given mode should show by default, for the given x-range. + void GetDefaultYRange(const PMusrViewMode &mode, Double_t xmin, Double_t xmax, + Double_t &ymin, Double_t &ymax) const; + /// Axis titles for the given mode. + void GetAxisTitles(const PMusrViewMode &mode, std::string &xTitle, std::string &yTitle) const; + + /// True if the msr-file PLOT block asks for a logarithmic x-axis. + Bool_t IsLogX() const; + /// True if the msr-file PLOT block asks for a logarithmic y-axis. + Bool_t IsLogY() const; + + //-- text blocks ------------------------------------------------------- + /// The msr-file title. + const std::string &GetTitle() const { return fTitle; } + /// The fit parameter block, one line per parameter. + const std::vector &GetParameterLines() const { return fParameterLines; } + /// The theory and functions block, one line each. + const std::vector &GetTheoryLines() const { return fTheoryLines; } + /// The run info/legend block, one line per data set. + const std::vector &GetInfoLines() const { return fInfoLines; } + /// The musrfit statistics line (date, chisq/maxLH, NDF). + const std::string &GetStatisticLine() const { return fStatisticLine; } + /// RRF information, empty if no RRF is used. + const std::string &GetRRFText() const { return fRRFText; } + + /// Ascii dump of what the given mode shows, restricted to [xmin, xmax]. + Bool_t ExportData(const char *fileName, const PMusrViewMode &mode, + Double_t xmin, Double_t xmax); + + //-- access for the front ends ---------------------------------------- + const PMusrViewDataList &GetData() const { return fData; } + const PMusrViewNonMusrDataList &GetNonMusrData() const { return fNonMusrData; } + const PMsrFourierStructure &GetFourier() const { return fFourier; } + const PIntVector &GetMarkerList() const { return fMarkerList; } + const PIntVector &GetColorList() const { return fColorList; } + + private: + Bool_t fValid{false}; ///< true if everything looks OK + std::string fLastError{}; ///< last error message + Bool_t fTheoAsData{false}; ///< calculate theory points only at the data points + Bool_t fScaleN0AndBkg{true}; ///< true=N0/background scaled to 1/ns, otherwise 1/bin + Int_t fPlotType{-1}; ///< MSR_PLOT_* tag + Int_t fPlotNumber{-1}; ///< plot block number + + Bool_t fFourierPrepared{false}; ///< Fourier transforms of data/theory are available + Bool_t fDifferencePrepared{false}; ///< time domain difference is available + Bool_t fDiffFourierPrepared{false}; ///< d-f, Fourier of the difference, is available + Bool_t fFourierDiffPrepared{false}; ///< f-d, difference of the Fourier spectra, is available + Bool_t fAveragePrepared{false}; ///< average is available + std::string fAverageKey{}; ///< mode key for which the average was calculated + + Double_t fXmin{0.0}, fXmax{0.0}; ///< time domain data range + Double_t fYmin{0.0}, fYmax{0.0}; ///< time domain value range + Bool_t fXRangePresent{false}; ///< x-range given in the msr-file + Bool_t fYRangePresent{false}; ///< y-range given in the msr-file + + PDoubleVector fCurrentFourierPhase{}; ///< currently applied Fourier phase(s) + + std::string fTitle{}; ///< msr-file title + std::vector fParameterLines{}; ///< fit parameter block + std::vector fTheoryLines{}; ///< theory/functions block + std::vector fInfoLines{}; ///< run info/legend block + std::string fStatisticLine{}; ///< musrfit statistics line + std::string fRRFText{}; ///< RRF information + + TH1F *fHistoFrame{nullptr}; ///< frame defining the visible time range + Bool_t fVisibleRangeSet{false}; ///< true once the user restricted the visible time range + +#ifndef __MAKECLING__ + PMsrHandler *fMsrHandler{nullptr}; ///< msr-file handler + PRunListCollection *fRunList{nullptr}; ///< data handler +#endif // __MAKECLING__ + + PMusrViewDataSet fDataAvg{}; ///< averaged data set + PMusrViewDataList fData{}; ///< all histogram data sets + PMusrViewNonMusrDataList fNonMusrData{}; ///< all non-muSR error graph data sets + + PMsrFourierStructure fFourier{}; ///< Fourier transform settings + PIntVector fMarkerList{}; ///< list of markers + PIntVector fColorList{}; ///< list of colors + + //-- data preparation, taken over unchanged from the former PMusrCanvas -- + void InitFourier(); + void InitAverage(); + void InitDataSet(PMusrViewDataSet &dataSet); + void InitDataSet(PMusrViewNonMusrDataSet &dataSet); + void CleanupDataSet(PMusrViewDataSet &dataSet); + void CleanupDataSet(PMusrViewNonMusrDataSet &dataSet); + void HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data); + void HandleNonMusrDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data); + void HandleDifference(); + void HandleFourier(const PMusrViewMode &mode); + void HandleDifferenceFourier(); + void HandleFourierDifference(); + void HandleAverage(); + void CleanupDifference(); + void CleanupFourier(); + void CleanupFourierDifference(); + void CleanupAverage(); + + void CalcPhaseOptReFT(); + Double_t CalculateDiff(const Double_t x, const Double_t y, TH1F *theo); + Double_t CalculateDiff(const Double_t x, const Double_t y, TGraphErrors *theo); + Int_t FindBin(const Double_t x, TGraphErrors *graph); + + Double_t GetMaximum(TH1F *histo, Double_t xmin = -1.0, Double_t xmax = -1.0); + Double_t GetMinimum(TH1F *histo, Double_t xmin = -1.0, Double_t xmax = -1.0); + Double_t GetMaximum(TGraphErrors *graph, Double_t xmin = -1.0, Double_t xmax = -1.0); + Double_t GetMinimum(TGraphErrors *graph, Double_t xmin = -1.0, Double_t xmax = -1.0); + + Bool_t IsScaleN0AndBkg(); + UInt_t GetNeededAccuracy(PMsrParamStructure param); + Double_t GetInterpolatedValue(TH1F *histo, Double_t xVal); + + void GetExportDataSet(const TH1F *data, const Double_t xmin, const Double_t xmax, + PMusrViewAsciiDumpVector &dumpData, const Bool_t hasError = true); + + //-- text block generation -------------------------------------------- + void BuildTextBlocks(); + void BuildRRFText(); + + //-- helper ------------------------------------------------------------ + /// Picks one histogram of a data set according to mode and role. + TH1F *SelectHisto(const PMusrViewDataSet &dataSet, const PMusrViewMode &mode, + Bool_t theory, Bool_t imagPart) const; +}; + +} // namespace musrview + +#endif // _PMUSRVIEWMODEL_H_ diff --git a/src/include/PMusrWebCanvas.h b/src/include/PMusrWebCanvas.h new file mode 100644 index 00000000..b5744054 --- /dev/null +++ b/src/include/PMusrWebCanvas.h @@ -0,0 +1,185 @@ +/*************************************************************************** + + PMusrWebCanvas.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2026 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 _PMUSRWEBCANVAS_H_ +#define _PMUSRWEBCANVAS_H_ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "PMusrViewMode.h" +#include "PMusrViewModel.h" + +namespace musrview { + +//-------------------------------------------------------------------------- +/** + * \brief ROOT7/web based front end of musrview. + * + *

Following the ROOT7 architecture (see S. Linev et al., J. Phys. Conf. Ser. + * 1525 (2020) 012092), the graphics are fully decoupled from the data: the + * C++ side is the server and holds the model (PMusrViewModel), the client is a + * web page which renders the spectra with JSROOT and the text blocks as plain + * HTML. Both talk to each other through an ROOT::RWebWindow websocket, the + * payload being JSON produced by TBufferJSON. + * + *

This buys three things which the former TCanvas based implementation could + * not provide: + * + * - The parameter/theory/info blocks are HTML panels with their own scrollbar + * instead of TPaveText/TLegend objects which scale their font down until the + * content is no longer readable. + * - Every block except the spectra can be shown and hidden individually. + * - The command keys are handled by the page, which knows whether a text input + * currently has the focus. Typing a 'd' or an 'f' into a text/latex label can + * therefore no longer trigger the difference/Fourier command. + * + *

The window can either be displayed in the user's standard browser or in an + * embedded Qt6 WebEngine window, in which case musrview looks and behaves like + * the stand-alone application it has always been. Which one is used is decided + * by the string handed to Show(). + * + *

The graphics export (png, pdf, ...) does not go through the browser but + * builds a plain TCanvas with the classic layout, so that batch operation stays + * independent of any browser installation. + */ +class PMusrWebCanvas +{ + public: + PMusrWebCanvas(Int_t plotNumber, const std::string &title, Int_t ww, Int_t wh, + PMsrFourierStructure fourierDefault, + const PIntVector &markerList, const PIntVector &colorList, + Bool_t batch, Bool_t fourier, Bool_t avg, Bool_t theoAsData); + virtual ~PMusrWebCanvas(); + + PMusrWebCanvas(const PMusrWebCanvas &) = delete; + PMusrWebCanvas &operator=(const PMusrWebCanvas &) = delete; + +#ifndef __MAKECLING__ + /// Hands the msr-file handler to the model. + void SetMsrHandler(PMsrHandler *msrHandler); + /// Hands the run collection to the model. + void SetRunListCollection(PRunListCollection *runList); +#endif // __MAKECLING__ + + /// Builds the data of this plot block. Has to be called before anything else. + Bool_t Build(); + /// True if everything is in a usable state. + Bool_t IsValid() const { return fValid; } + /// Last error message. + const std::string &GetLastError() const { return fLastError; } + + /// Seconds after which the window closes by itself; <= 0 means no timeout. + void SetTimeout(Int_t timeout) { fTimeout = timeout; } + + /// Opens the window. \p where is "" (standard browser), "qt6", "chrome", ... + Bool_t Show(const std::string &where); + /// True once the user has closed the window. + Bool_t IsDone() const { return fDone; } + /// Closes the window and drops the connection. Needed so that musrview can + /// terminate: an open RWebWindow keeps the process alive. + void Close(); + /// Called when the user quits; used to terminate the application. + void SetDoneCallback(std::function cb) { fDoneCallback = std::move(cb); } + + /// Writes the classic canvas layout into a graphics file. + Bool_t SaveGraphics(const std::string &msrFileName, const std::string &graphicsFormat); + /// Writes the currently shown data as ascii. + Bool_t ExportData(const std::string &fileName); + + /// The current view mode. + const PMusrViewMode &GetMode() const { return fMode; } + /// The underlying model. + PMusrViewModel *GetModel() { return fModel.get(); } + + private: + Int_t fPlotNumber{-1}; ///< plot block number + std::string fTitle{}; ///< msr-file title + Int_t fWidth{800}; ///< window width + Int_t fHeight{600}; ///< window height + Bool_t fBatchMode{false}; ///< no window is opened + Bool_t fStartWithFourier{false}; ///< start with the Fourier representation + Bool_t fStartWithAvg{false}; ///< start with the averaged representation + Bool_t fValid{false}; ///< true if usable + Bool_t fDone{false}; ///< true once the user quit + Int_t fTimeout{0}; ///< timeout in seconds + std::string fLastError{}; ///< last error message + + std::unique_ptr fModel; ///< the data layer + PMusrViewMode fMode{}; ///< what is currently shown + + std::shared_ptr fWindow; ///< the web window + std::function fDoneCallback{}; ///< invoked when the user quits + + std::unique_ptr fFrame; ///< drawing frame of the spectra + std::unique_ptr fMultiGraph; ///< non-muSR data/theory + std::unique_ptr fPlotLegend; ///< legend of the non-muSR plot + std::unique_ptr fPhaseLatex; ///< Fourier phase label + std::unique_ptr fRRFLatex; ///< RRF label + + //-- web communication ------------------------------------------------- + void ProcessMessage(unsigned connid, const std::string &msg); + void SendState(unsigned connid); + void SendPlot(unsigned connid); + void SendAll(unsigned connid); + std::string BuildStateJson() const; + /// Describes each run by its info line and its legend symbol. + std::string BuildRunsJson() const; + + //-- plotting ---------------------------------------------------------- + /// Collects the objects of the current view, each with its draw option. + Bool_t CollectDrawables(std::vector> &out); + /// Draws the spectra of the current mode into the given pad (graphics export). + Bool_t FillPlotPad(TVirtualPad *pad); + /// Serializes the objects of the current view for JSROOT. + Bool_t BuildPlotJson(TString &json); + + //-- helper ------------------------------------------------------------ + static std::string JsonEscape(const std::string &str); + static std::string JsonArray(const std::vector &vec); + /// Applies a command coming from the page. Returns true if a redraw is due. + Bool_t ApplyCommand(const std::string &cmd); +}; + +} // namespace musrview + +#endif // _PMUSRWEBCANVAS_H_ diff --git a/src/musrview.cpp b/src/musrview.cpp index dd8a1ac1..2e3e7178 100644 --- a/src/musrview.cpp +++ b/src/musrview.cpp @@ -31,33 +31,35 @@ #include "config.h" #endif -#include -#include -#include +#include +#include +#include +#include #include #include +#include #include #include -#include - #include -#include #include +#include +#include + +#include #ifdef HAVE_GIT_REV_H #include "git-revision.h" #endif #include "PMusr.h" -#include "PStartupHandler.h" +#include "PMsgBox.h" #include "PMsrHandler.h" +#include "PMusrWebCanvas.h" #include "PRunDataHandler.h" #include "PRunListCollection.h" -#include "PMusrCanvas.h" -#include "PMsgBox.h" - +#include "PStartupHandler.h" //-------------------------------------------------------------------------- /** @@ -74,8 +76,16 @@ void musrview_syntax() std::cout << std::endl << " -f, --fourier: will directly present the Fourier transform of the ."; std::cout << std::endl << " -a, --avg: will directly present the averaged data/Fourier of the ."; std::cout << std::endl << " -1, --one_to_one: calculate theory points only at data points."; + std::cout << std::endl << " display options (musrview is web based, see the musrfit documentation):"; + std::cout << std::endl << " --embedded : show the plots in an embedded browser window, i.e. musrview"; + std::cout << std::endl << " behaves like a stand-alone application. This is the default"; + std::cout << std::endl << " as long as ROOT provides the Qt6 web display."; + std::cout << std::endl << " --browser : show the plots in the standard web browser of the system."; + std::cout << std::endl << " --browser= : use a specific display, e.g. chrome, firefox, qt6, server."; + std::cout << std::endl << " 'server' only prints the URL and starts no browser at all,"; + std::cout << std::endl << " which is what to use when musrfit runs on a remote machine."; std::cout << std::endl << " --: "; - std::cout << std::endl << " will produce a graphics-output-file without starting a root session."; + std::cout << std::endl << " will produce a graphics-output-file without opening a window."; std::cout << std::endl << " the name is based on the , e.g. 3310.msr -> 3310_0.png"; std::cout << std::endl << " supported graphic-format-extension:"; std::cout << std::endl << " eps, pdf, gif, jpg, png, svg, xpm, root"; @@ -91,24 +101,64 @@ void musrview_syntax() } //-------------------------------------------------------------------------- +/** + *

Shows an error message in a message box. + * + * \param errMsg the message to be shown + */ void musrview_error_msg(std::string errMsg) { - int argc=0; - char **argv; + int argc=0; + char **argv=nullptr; TApplication app("musrviewErrorMsg", &argc, argv); new PMsgBox(errMsg, gClient->GetRoot(), 600, 200); app.Run(); } +//-------------------------------------------------------------------------- +/** + *

Decides which web display musrview should use. + * + *

Unless the user asked for something specific, the embedded Qt6 WebEngine + * window is preferred, since it makes musrview appear as the stand-alone + * application it has always been. If ROOT was built without the Qt6 web + * display, the standard browser of the system is used instead. + * + * \param requested what the user asked for, empty if nothing was requested + * + * \return the display string to be handed to ROOT::RWebWindow::Show() + */ +std::string musrview_select_display(const std::string &requested) +{ + if (!requested.empty()) + return requested; + + if (const char *env = gSystem->Getenv("MUSRVIEW_BROWSER")) + if (strlen(env) > 0) + return std::string(env); + + // prefer the embedded window if ROOT provides it + gSystem->RedirectOutput("/dev/null", "a"); + const Int_t loaded = gSystem->Load("libROOTQt6WebDisplay"); + gSystem->RedirectOutput(nullptr); + + if (loaded >= 0) + return "qt6"; + + std::cout << std::endl << ">> musrview **INFO** ROOT has no Qt6 web display, using the standard browser."; + std::cout << std::endl << ">> Rebuild ROOT with -Dqt6web=ON for an embedded window." << std::endl; + + return std::string(); +} + //-------------------------------------------------------------------------- /** *

The musrview program is used to show muSR fit results in graphical form. * From it also Fourier transforms, difference between data an theory can be formed. - * For a detailed description/usage of the program, please see - * \htmlonly musrview online help - * \endhtmlonly - * \latexonly musrview online help: \texttt{http://lmu.web.psi.ch/musrfit/user/html/user-manual.html#musrview} - * \endlatexonly + * + *

Since version 1.13 musrview is based on the ROOT7 web graphics: the plots + * are rendered by JSROOT either in an embedded browser window or in the + * standard browser of the system, while the C++ side only holds the data. * * return: * - PMUSR_SUCCESS if everthing went smooth @@ -131,6 +181,7 @@ int main(int argc, char *argv[]) bool graphicsOutput{false}; bool asciiOutput{false}; std::string graphicsExtension{}; + std::string display{}; int timeout{0}; bool show_errMsgBox{false}; std::stringstream errMsg; @@ -182,6 +233,16 @@ int main(int argc, char *argv[]) avg = true; } else if (!strcmp(argv[i], "-1") || !strcmp(argv[i], "--one_to_one")) { theoAtData = true; + } else if (!strcmp(argv[i], "--embedded")) { + display = "qt6"; + } else if (!strcmp(argv[i], "--browser")) { + display = "default"; + } else if (!strncmp(argv[i], "--browser=", 10)) { + display = argv[i]+10; + if (display.empty()) { + show_syntax = true; + break; + } } else if (!strcmp(argv[i], "--eps") || !strcmp(argv[i], "--pdf") || !strcmp(argv[i], "--gif") || !strcmp(argv[i], "--jpg") || !strcmp(argv[i], "--png") || !strcmp(argv[i], "--svg") || !strcmp(argv[i], "--xpm") || !strcmp(argv[i], "--root")) { @@ -217,6 +278,8 @@ int main(int argc, char *argv[]) return PMUSR_WRONG_STARTUP_SYNTAX; } + const bool batchMode = graphicsOutput || asciiOutput; + // read startup file std::string startup_path_name{}; std::unique_ptr saxParser = std::make_unique(); @@ -227,11 +290,9 @@ int main(int argc, char *argv[]) } else { startup_path_name = startupHandler->GetStartupFilePath().Data(); saxParser->ConnectToHandler("PStartupHandler", startupHandler.get()); - //status = saxParser->ParseFile(startup_path_name); - // parsing the file as above seems to lead to problems in certain environments; - // use the parseXmlFile function instead (see PStartupHandler.cpp for the definition) + // parsing the file via ParseFile() seems to lead to problems in certain + // environments; use parseXmlFile (see PStartupHandler.cpp) instead status = parseXmlFile(saxParser.get(), startup_path_name.c_str()); - // check for parse errors if (status) { // error std::cerr << std::endl << ">> musrview **WARNING** Reading/parsing musrfit_startup.xml failed."; std::cerr << std::endl << ">> Any graph will appear with random symbols and colors!"; @@ -265,13 +326,13 @@ int main(int argc, char *argv[]) musrview_error_msg(errMsg.str()); return status; } + // make a plot list vector PMsrPlotList *msrPlotList = msrHandler->GetMsrPlotList(); PIntVector plotList; bool runPresent; for (unsigned int i=0; isize(); i++) { for (unsigned int j=0; jat(i).fRuns.size(); j++) { - // check that run is not already in the plotList runPresent = false; for (unsigned int k=0; kat(i).fRuns[j] == static_cast(plotList[k])) { @@ -279,9 +340,8 @@ int main(int argc, char *argv[]) break; } } - if (!runPresent) { + if (!runPresent) plotList.push_back(static_cast(msrPlotList->at(i).fRuns[j])); - } } } @@ -306,10 +366,8 @@ int main(int argc, char *argv[]) // generate the necessary histogramms for the view std::unique_ptr runListCollection; if (result == PMUSR_SUCCESS) { - // feed all the necessary histogramms for the view runListCollection = std::make_unique(msrHandler.get(), dataHandler.get(), theoAtData); for (unsigned int i=0; iGetMsrRunList()->size(); i++) { - // if run is in plotList add it, otherwise go to the next runPresent = false; for (unsigned int j=0; j(plotList[j]) == i+1) { @@ -332,110 +390,140 @@ int main(int argc, char *argv[]) } } - if (result == PMUSR_SUCCESS) { - // create the ROOT application needed for PMusrCanvas - // and pass it only essential arguments - std::vector args; - args.push_back(argv[0]); // program name - if (graphicsOutput || asciiOutput) { - args.push_back((char*)"-b"); // batch mode flag - } - int cc = args.size(); - TApplication app("App", &cc, args.data()); - - std::vector canvasVector; - PMusrCanvas *musrCanvas = nullptr; - - bool ok = true; - for (unsigned int i=0; iGetMsrPlotList()->size(); i++) { - - if (startupHandler) - musrCanvas = new PMusrCanvas(i, msrHandler->GetMsrTitle()->Data(), - 10+i*100, 10+i*100, 800, 600, - startupHandler->GetFourierDefaults(), - startupHandler->GetMarkerList(), - startupHandler->GetColorList(), - graphicsOutput||asciiOutput, - fourier, avg, theoAtData); - else - musrCanvas = new PMusrCanvas(i, msrHandler->GetMsrTitle()->Data(), - 10+i*100, 10+i*100, 800, 600, - graphicsOutput||asciiOutput, - fourier, avg, theoAtData); - - if (musrCanvas == nullptr) { - errMsg << ">> musrview **SEVERE ERROR** Couldn't invoke all necessary objects, will quit.\n"; - std::cerr << errMsg.str(); - if (show_errMsgBox) - musrview_error_msg(errMsg.str()); - ok = false; - result = PMUSR_MSR_RUN_ERROR; - break; - } - if (!musrCanvas->IsValid()) { - errMsg << ">> musrview **SEVERE ERROR** Couldn't invoke all necessary objects, will quit.\n"; - std::cerr << errMsg.str(); - if (show_errMsgBox) - musrview_error_msg(errMsg.str()); - ok = false; - result = PMUSR_MSR_RUN_ERROR; - delete musrCanvas; - break; - } - - musrCanvas->SetTimeout(timeout); - - // ugly but rootcling cannot handle the spirit-parser framework - musrCanvas->SetMsrHandler(msrHandler.get()); - musrCanvas->SetRunListCollection(runListCollection.get()); - - musrCanvas->UpdateParamTheoryPad(); - musrCanvas->UpdateDataTheoryPad(); - musrCanvas->UpdateInfoPad(); - - if (!musrCanvas->IsValid()) { // something went wrong - ok = false; - result = PMUSR_MSR_RUN_ERROR; - break; - } - - musrCanvas->Connect("Done(Int_t)", "TApplication", &app, "Terminate(Int_t)"); - - if (graphicsOutput) { - musrCanvas->SaveGraphicsAndQuit(fileName.c_str(), graphicsExtension.c_str()); - } - - if (asciiOutput) { - // generate export data file name - TString str(fileName.c_str()); - str.Remove(str.Last('.')); - str += ".dat"; - // save data in batch mode - musrCanvas->ExportData(str.Data()); - musrCanvas->Done(0); - } - - // keep musrCanvas objects - canvasVector.push_back(musrCanvas); - } - - // check that everything is ok - if (ok) - app.Run(true); // true needed that Run will return after quit so that cleanup works - - // clean up - char canvasName[32]; - for (unsigned int i=0; iGetListOfCanvases()->FindObject(canvasName) != nullptr) { - canvasVector[i]->~PMusrCanvas(); - } - } - canvasVector.clear(); + if (result != PMUSR_SUCCESS) { + plotList.clear(); + return result; } - // clean up + // A ROOT application is needed so that the web display and the graphics + // export find a working environment. + // + // In the interactive case no TCanvas is created at all - the spectra travel + // to JSROOT as a list of drawables - hence ROOT never opens a window of its + // own. Only the graphics export builds a TCanvas, and that runs in batch. + std::vector args; + args.push_back(argv[0]); + if (batchMode) + args.push_back(const_cast("-b")); + int cc = static_cast(args.size()); + TApplication app("musrview", &cc, args.data()); + + std::vector> canvasList; + + bool ok = true; + for (unsigned int i=0; iGetMsrPlotList()->size(); i++) { + auto canvas = std::make_unique( + static_cast(i), msrHandler->GetMsrTitle()->Data(), + 1200, 900, + startupHandler ? startupHandler->GetFourierDefaults() : PMsrFourierStructure{}, + startupHandler ? startupHandler->GetMarkerList() : PIntVector{}, + startupHandler ? startupHandler->GetColorList() : PIntVector{}, + batchMode, fourier, avg, theoAtData); + + // ugly but rootcling cannot handle the spirit-parser framework + canvas->SetMsrHandler(msrHandler.get()); + canvas->SetRunListCollection(runListCollection.get()); + + if (!canvas->Build()) { + errMsg << ">> musrview **ERROR** " << canvas->GetLastError() << "\n"; + std::cerr << errMsg.str(); + if (show_errMsgBox) + musrview_error_msg(errMsg.str()); + ok = false; + result = PMUSR_MSR_RUN_ERROR; + break; + } + + canvas->SetTimeout(timeout); + canvasList.push_back(std::move(canvas)); + } + + if (!ok) { + canvasList.clear(); + plotList.clear(); + return result; + } + + if (graphicsOutput) { + for (auto &canvas : canvasList) { + if (!canvas->SaveGraphics(fileName, graphicsExtension)) { + std::cerr << std::endl << ">> musrview **ERROR** " << canvas->GetLastError() << std::endl; + result = PMUSR_MSR_RUN_ERROR; + } + } + } + + if (asciiOutput) { + // generate the export data file name + TString str(fileName.c_str()); + str.Remove(str.Last('.')); + for (size_t i=0; i 1) { + name += "_"; + name += static_cast(i); + } + name += ".dat"; + if (!canvasList[i]->ExportData(name.Data())) { + std::cerr << std::endl << ">> musrview **ERROR** " << canvasList[i]->GetLastError() << std::endl; + result = PMUSR_MSR_RUN_ERROR; + } + } + } + + if (!batchMode) { + const std::string where = musrview_select_display(display); + + // one window per PLOT block, exactly as before + ROOT::RWebWindowsManager::SetSingleConnMode(false); + + for (auto &canvas : canvasList) { + if (!canvas->Show(where)) { + std::cerr << std::endl << ">> musrview **ERROR** " << canvas->GetLastError() << std::endl; + result = PMUSR_MSR_RUN_ERROR; + ok = false; + break; + } + } + + if (ok) { + // run until every window has been closed, or until the timeout expired + const auto start = std::chrono::steady_clock::now(); + while (true) { + gSystem->ProcessEvents(); + gSystem->Sleep(20); + + bool allDone = true; + for (auto &canvas : canvasList) { + if (!canvas->IsDone()) { + allDone = false; + break; + } + } + if (allDone) + break; + + if (timeout > 0) { + const auto now = std::chrono::steady_clock::now(); + if (std::chrono::duration_cast(now-start).count() >= timeout) { + std::cout << std::endl << ">> musrview: timeout of " << timeout << " s reached, terminating." << std::endl; + break; + } + } + } + } + + // End the process here. Tearing down an RWebWindow which drives the + // embedded display blocks, because the Qt event loop ROOT started on our + // behalf is never stopped again. Everything of interest has been written by + // now, so leave immediately rather than hang after the last window closed. + std::cout.flush(); + std::cerr.flush(); + gSystem->Exit(result, kFALSE); + } + + // only reached in batch mode, where no window was ever opened + canvasList.clear(); plotList.clear(); return result; diff --git a/src/musrview_web/musrview.html b/src/musrview_web/musrview.html new file mode 100644 index 00000000..fc7d9292 --- /dev/null +++ b/src/musrview_web/musrview.html @@ -0,0 +1,705 @@ + + + + + musrview + + + + + +

+ +
+
+ + +
+ +
+ + + + +
+ + + + +
+ + + + +
+ +
+ + + + + +
+ +
+ + + +
+ +
 
+ +
+ +
+
+

Fit parameters

+ 

+      

+      

+        
Theory

+        

+      

+    

+
+    

+      
Run info

+      

+    

+
+    
connecting …

+
+  

+
+
+
+
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 71813837..6a31c407 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
 #--- musrfit integration tests ------------------------------------------------
 add_subdirectory(maxLH_check)
+add_subdirectory(viewmode_check)
 add_subdirectory(musrview_check)
diff --git a/tests/musrview_check/CMakeLists.txt b/tests/musrview_check/CMakeLists.txt
index 387f1729..e2a2e4a2 100644
--- a/tests/musrview_check/CMakeLists.txt
+++ b/tests/musrview_check/CMakeLists.txt
@@ -60,6 +60,25 @@ macro(add_musrview_test test_name msr_file)
   )
 endmacro()
 
+# add_musrview_test_dir(sub_dir test_name msr_file [extra_musrview_args...])
+# Like add_musrview_test(), but for msr-files whose data lives in a
+# sub-directory of doc/examples; that sub-directory also becomes the working
+# directory so that the relative data paths of the msr-file resolve.
+macro(add_musrview_test_dir sub_dir test_name msr_file)
+  add_test(
+    NAME ${test_name}
+    COMMAND ${MUSRVIEW_CHECK_PYTHON3}
+      ${CMAKE_CURRENT_SOURCE_DIR}/musrview_check.py
+      $
+      ${CMAKE_SOURCE_DIR}/doc/examples/${sub_dir}/${msr_file}
+      ${CMAKE_CURRENT_SOURCE_DIR}/ref
+      ${test_name}
+      --tol 0.03
+      ${ARGN}
+    WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/doc/examples/${sub_dir}
+  )
+endmacro()
+
 #--- doc/examples tests -------------------------------------------------------
 
 # single histogram fits
@@ -170,3 +189,13 @@ add_musrview_test(musrview-histo-PSI-BIN-scale-ufr-lc-fourier-avg
 add_musrview_test(musrview-histo-PSI-BIN-scale-ufr-no-lc                           ViewOpts/test-histo-PSI-BIN-scale-use_fit_ranges-no-lifetimecorrection.msr)
 add_musrview_test(musrview-histo-PSI-BIN-scale-ufr-no-lc-fourier                   ViewOpts/test-histo-PSI-BIN-scale-use_fit_ranges-no-lifetimecorrection.msr -f)
 add_musrview_test(musrview-histo-PSI-BIN-scale-ufr-no-lc-fourier-avg               ViewOpts/test-histo-PSI-BIN-scale-use_fit_ranges-no-lifetimecorrection.msr -f -a)
+
+#--- non-muSR plot ------------------------------------------------------------
+#--- Covers the only fit type whose axis labels come from the data file rather
+#--- than from the plot type. Needs the BMW user function library, hence the
+#--- guard, and the freshly built library has to be findable at run time.
+if (BMWlibs)
+  add_musrview_test_dir(BMWlibs musrview-nonMusr-GapIntegrals test-libGapIntegrals-ASCII.msr)
+  set_tests_properties(musrview-nonMusr-GapIntegrals PROPERTIES
+    ENVIRONMENT "LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}/src/external/libGapIntegrals:$ENV{LD_LIBRARY_PATH}")
+endif (BMWlibs)
diff --git a/tests/musrview_check/ref/musrview-nonMusr-GapIntegrals/test-libGapIntegrals-ASCII_0.png b/tests/musrview_check/ref/musrview-nonMusr-GapIntegrals/test-libGapIntegrals-ASCII_0.png
new file mode 100644
index 0000000000000000000000000000000000000000..6c7ea008ab03f80118e06b2666509ba49555873e
GIT binary patch
literal 17383
zcmb_^bwHF&_x244sFXf}G#KPci7X(ENGuKR(j{RK0t!ll1qvb!N=xt3B}yX#A|kn@
zw8YZ22uppl#`8XJ{Qmp)54kt)nYm}qoO4~*IpcF(ZB?4%EXN@TqCu!#)rTN5HV7h>
zIZ6sX*{qcs0e>j0HB_&Hzu>i_ROJsrXCTB?MWjdK(n!a>VU_r8|4GyLA5-o|R1OQ}
z+6#G?-b?xj@n-L++B&(JzRwcC=0P))i;Akl<=tb&2}6VB&=+LEDJLvP)|M6euTg#u
zajA80k4SBuDN4l>MYqQ*{fgX=&a6To!-7Uamm7%vvQkkzZm=$%lcEN1EXfs%j06tq
z3L3$iIkHF?_(P?J0sjij37~`1_->`(M
z7!}FsM7`xga4m8nDqoUAR%5*8`2%VW>oSIiLDvVihVa-1OXFD2XwKajQc}nOYpid#
zXJiQHcV*j<=71oNmV|kE{wQ51WA!%P7oIP?V!05TaeDOFIJTV6EafUT2i<}2CKyLe
zfx7d{9@|sF8*>*>8o4_jIPudG?XMc18C%|PeMq22O{-+?Vld)w{Vt^O>63>>wc2KO
zoHm4A>YrUBEMs2x4CFS~VvbdLeDgx?Q9#?$RFpYZlafNUtSzp$zu%4*O5&wfvSGk4
z{Dne=6I9h7q!Lo&uW9j{--y!cob5V)ku-+)t)9~IAIn0q60?`yZy<{du90sdHg#vH`%7qD+m|y;
zj$bj8=}JNFSTE)x$gV%9}|%ErL+J?vg~dboCr{9-x1934#WO}v_|byT}p%2c7W^()Vx-l!TZMOMAx3jG;8nqJH@jz=Zbwf1kMh&XT^=$9
z%&6OKf=~dXk2=O4CaN6C)qVHVr-D8qc|LLHF|GuDxnax7`pks)-6n3R2+8BsRqnFRD{p{3%CUi(Z*66odX8MNpS2R~*aVpKUCxt6nN
zUjxTzKSokz$E)tsz$M(L1k;!ET5t5}a!L>;yS)wdjUPqWVZ!oERP-ymXWCx{vtn&8
zsK@l(y)#gsYyTR4xzc0=^FDoa^kv3zNahB5dRGN0(v*Gpf^z!F#k#dKm0zy;Jjvhw
zVve@%pvt7HT-Yrdh^)kDklTFi<~@}wUOInKYHI2_6H9`n&+)?Aif0(w#Nd1L7vQm3
zh5ZAW&9$(RV{lVFLb%l$lZb)+m%i_Vm4qS<=~BMa(S1*xe413g$m~->wX0bgPXm*r
zV!!H|-N3q;`<{%8&#zs?z|f-6_d_~uFD0yS`(APo&#%1@gemV|dD1P!H7`-rYJhYv
z#M
zC=(OfVz)~R&J|7HGRqVS1}d^Inr>uZ{CIme84DTrn5xoR4WdsM!3Q5YY?9jfIAo6a
zGSm$cw}`Ln!D
zv3HC77OqDEjWZtpbf=UNND>G#Qqadr`Bv8z*nl@qYlUy5ASmFPk5djU1dX5h<#$Qv
z4O+bn_C*;*3v+=`7py0Qgc$F5U8R>W1Ia>K)-RbsMmRwChbTCuTN
zx7k|oz)D|bp}o*=`Sz8)`GfKLIN|N}C-U-|m0PnU(5fMpz08wv
zZc@@(w4>S5nRiWA|WA`pU+faz1>4EG@5$g)y%Yf
zgrdQa`Mh9^du#9ut~{8-x5SKVu)_;qQ^k!%dCP}Gysb7eg6HKSXz^K^+MS80)~lk~
z^7}ZQ^&cJn2c=&&ED!OPcj@G$r8Oe2E-%1LrKE{3s>aoY!a5-#OZ#_Rg5j`y$(Nb)
z>exq3u53pdJY>XY422CbWAw(cD3qD{afTy`$+PRC)Qtjd>oO7s88>cqbQ7dsUzg!n
zz~yPuJ*t28j)J&+i{TbS3ATcf=}ex@SS&Wq(vLatx)z0vTBNP$Yc-)v7CJ<%!BDeW
zxh#xF7CxU4&y6^DJMBxYZzLgWwD!X7t_ryTC2gzkidm>{S{^$soD3$Jv~8cx8+44c2Gm?zgt`SLR1!vs->xZIpD-IR0E>RC*d4PF;N
z4KeT1K}8|9W)TJH1IgyYqk=ZCSRm;BBcvnNO$xg3$&>`DdKSrQboHkR!igpbo&GP|
zl)uzWVFBwYV>5{>5y&0mB?NTZLpR`dnP(IGshqJSG59$9&RbMEs(1P(PD79dY?K6>
zQ%}>P*5^I+En`gWhoBbIW&zypnJjicYj3JI;B)uQns{U1`%%u}#fM-fVyIV3x{bI+
zyHn;T%>u(u)!^u1h`-MJisdNY?cLL-%?Qkj;1_(#JN1|P56#+}K2pCaPnIl70m&k
zGdHg4_4=?as5=;Pl0xH!lkB>MjrNj+Lc+7MhHllXcO?%B&I`tHEnqER8u
z3Ee`2^_T}(iuvJ8|EuDI=;FignWJ$rK93W^_Eo_Q3@D?1c$}xPe|KU>)%^hd;5_=y
zR;Kx5mJL=aw~qkfjGb$Be(_Sw&YHCt&8?A)OjNYXQ5njJgbg0KfMCD^$UAV6BW!v!
zan`d~T)lZSM!LewQ2(J!$1n$DA-AZ@sCBU}{LR#wB;qfd5EJ%m#<%6REF-z>Gvwi}
zG91oK5OEUfA_d-fmSTX%o0Kg0e0Gzyu56M2D>^J`R*s?0(FAAoFw;0Izr0XtD@8Nn
zs;-6lm(zhh7F8Kk5aiFW65xe5i;6$vo2T%4>{XhuZsBeHAC|s;>7=ev7`_noTmgKa
z8i6NAr}b)$XhCMBeEF(1gt9qgsEh9A&I=~mUsXu)q18#7L!x6>0+aG@F8P0~O`VC9
zpNp*Y)Ouw`D=i}O>dJ5rQ}MecU*cK|4&^lpLQpl*o$}lFnVdlUQqNMkpMu?+2N}1i
ziv{*wq~NJ{H7<7Q{Cw_Wb%S{M?Wgpw)!saNN}`ik?51R^sf52P!_hXgFm0ND`r6a$
zd>{-ufL>H+f~Vds2|op4J|IPlCn_c|Cqa`(U-vB9$oVhyOo>3-MzwvXEL?o-zRD#>
zxzxL+Ll=%a=u~+hfr=^%R;VvrVBZy{R(k(Kg4wIkD4o8gnA&4@iWaR`dRf0Jo0H6b
zrIj9aq2C^#QoBL3UOpdD#`t(&P_w6(>8w7+L_@+|R#w`Gr2c`NMQpBhEd84aSq))d
zy$~!b)6N}Sz~?JtWb@}m>%0a><&1ADOvesz>RxUuS@RHZRn@aaW~m>06JvivrQrya
zz3pMVb@7Iz25q*EqaZ+Pt+tU_x3!uu_C1DhRuBoZtH1Q@c#=TcYilKoHeXy~*|doW
z1G0M9o}td7jT06W0!*yluJ5>fh_6q%VF}DPk3lb-3QAfnhPAi6nh^8PZAsps$>Kc3P89e6G
zkUMkAjxM84mnwo8SvbCaH;N=|gjx4}6GPbuvkX*!k~ndZFlKbt#1*9FXj6odoX@tE
z@W4w!zz(dGYkKeTz#rEd;vX^j2BgBRucddY@fflsB-i8@2lLo&nBw#mgU9DaPeJ2x
zVp@h663xh~-Q+dPxcPb5%PFwoUI4l=aCVvx2cs%_7*Qq-q4-(wP!@r+iXp+hdFL*Hbp<(Q$lRRmo!1&$vyP
zXr#oiFD~XK)mPYqrPhzgX6#|w$Ys$-A#N_(;_ohrm8FouoN1Cmh$yI9V3+AqqTTxx
zm`zt^k03{FX72bV2rTBsDM*A(jpR#HryT8rBM_8mUAKDeF4}-~;MxtUe=heH9Y0Rb
zccV$&#K7d6-GZjb{E5p|kd>?)51(8!_094aBw90P;JVMoaYfmYV~w4BKT5oLTr5cH
z1AVH)37~VVwSM_^B|Km7x{OW_Qe~bwti_dWfu99<8_A{?Ds$_A%dBVrhOaGRT(ZqClmB20@*OG$}|n%UPD9Ds9>PK&v2|3T#2YL~oh
zh<&=vF3ydxKHWCxc;6kcllE@0q>foiO3LpbidxAEi6i5<)$Z~z#}t>fw>j$Q=-ezF
z+}YWAOi5?^d4GR@ZEej((r06)okxQrfc?7-DxNqfC-(z>hWm;q=FIy;tb6c0$R8V(
zm>Z_bnX_MZf2~*oKjsT2M!tHMAs8Cwn>$g-b@m$bt1{Dm`KpPI
z5bG;E8VMb<>0P(P6GBze%Nc1Q=}#^CoSZg+r9<{P9Ljp_E{`=34zLukR3!~HFc!q;$R%lW%KOQ%(avt{xJb%?mz0#(-gfcmq5aNj&8UORhbw8s1JpyoK<|6taiL)MFH1#?ex
z{)nOai>w}8Yilc*#Fgb`mGrKOi3v9;lm-O^$@ZxyTR4)(a6NCz(o5K@4!FGP?h-UN`
zT5__pqxEwz+^TwR4<7Vel|uD~e{B5~k|AXsIQrA3PKCZaNcRZnN^me|Br4(+
z<>Q7h`Gf7buJlZetlbNc?9A`U38BGqKZQlsc2PclIC`BRju=a>sHo^jmhH)wz-?`e
zR6Wp1meDKs4r#h#|ARx`PvBBsjV%Oqk_~_DF`O^LeCClC*4EaZ*#IWHlbzxRySBWt
zLPtk;_{0D3hqM1~f|tEcVgD(7xi>~-9em#%w>OW?1xVFzqg`S_nF_L4E$(dr&N%OF
zAR5i6tjBAI{M^;mwY0PZEFX`@`^eSR*2>Aq>}?G@`)?$SnVFln-zXK+0?Xde&;ZsZ
zCB9?!!Gj08ySvTZTR(nG=QmleO@W6!_)+J-I{t98$;H=qAM_cgPQ)$O!9AK)zy%MZ{s;3cLNIwYIU
zURi4k@3`-GtP95%K)Yi!KShi?(~&Sl_Moa4mM7B43bm(YlGal$F2T}oBF(gh{=maUAJNIG+iZ^X3LR&R+y6={l
z*FK@#;tMdmxAlI>$y|urN_Y`!C+z+p!Sm#ihhk~=>%Ga;jN6;IYn}w1K&rj%wW+42
zr}nuP*Q|SsD7)aDV2EFT57A&@ABi*f9TJQ$d4&`_9m0e6`j8`mY$PuzQ`u>myp!Fd
z2PJJ>R9pk`)?qlw-=qvNo6>yubcix;BLJ_n;uSpyL^!7`eOifJ+qv`g#U2|UUj~tMeqygq~
zv4m~lvdh83X&v$J_$)e&s2OxL7y>^bH*8wg#JSN!0ZSXAI?hwG%36j;&1KQqK}gQV
znJ{kBdLADEJn`Y$ZYa6qWEJrV;daKkBi{E-5f#tz$-NF1f}h}|3cNnm+r&QX1kv3a
zysZU8a{9v-ei|dU?;K<#+ZtvLe-(j4QosF5Y`HkfGx?@54=#FiBD)Salpo()qc%1+
zW;;^CRMKyDBqyb$$ompDmxe161nzI};_Xe$@_!gvvo19QIIMMg-A#7z@cnk~ek7(F
z_ts#@;z9iBlwOC?NsPJ%DiOLL{@ZgHUqA;>NggHTPgn18lP2c~5jBXFN$CLo37C5*
z4jBG=9}iW5RzBz`1Tf;#BJbTyjoqNw@#E(J&4R`fXL(%H-Mmwlyc;I>@
zoJZAOw}TBKG3Hac|H57>DLFZLd1(pZGy+^MOa7ob;P+a}Se2sW;z>rIv(YA7O4~
zN3j&dN}WNGYk+Nk09;3Z&SZu#*i7o>+lxw(T*k4}o<0J)%>ke2`oO-qZx3skAbNoBh|sj{ZrpK1GGFg*H@kH5B
zgO{+=w%aI6{0r~Waz{D{S~1eV0C1qmqf&Pa$N=^+vzKd_^D~xc)f{@W4IFNIyE&0V
za}`TiK2x@a8>+Am!>rsc+$uI(I5#g4hBk4;I0_&Gh|6tlZQ#Mh7unoZYo)FFw>|8M
zSPM&x(ssrML+LT2GrkYLJ~2=|W+oijWL>QCRw7Az44Cl2{%%!ORo5*pPS4G`uD2g3
zJB2uf+lT?|ay*iO@`Pq2SI0eB?XOP?^7>ytb#uahZ#~8TtrB?1ckjaMzM~Z6f5Y~+
zk&MqLjy&vwlij>B%T;B+E$G-?9wTk%4%z?tMOcRVUHb
zIfK}p@p3kXraSBwPVI}Ljroz<3$gbbimvw?P>UXgR-g4CupTlmr`~PmE8CccdCd~%
z=Ydnxna6-7zBAqC4O`Y+?M{}x9%#)ohWmni*7CKy@U_K%B$)G)Qa$VCg)!(Yze5RsZZQJ8^k&xyEh@kNcghMzGcMb_{LkGcVMaL3iWd`Xy62b!K~}LeN(t36wvNz+@naY
zRR{OK(s?@YBZ>I;-;>cFghn^KE;=78Zd38}tZM+x?-7OwF^pFRCJ8&%
zHaFGStd!*5>Z}tt*@`$SMAV$}snadFPOW=C2bmb>DJ3s_W>JpglL)+i=r@&^{viTa
zXy=ny$W%46Ab}EA%L#P4FNW&D_PBfZoX!$}w*W)+T$F;WK(?!s(o5+AMIs60_5pMZ
z`=B?yxrwBR)l)r6_7yT6|NW$aXhk|S$<56`18PaKvi9GWB8F5Qn&H%uV+gR6@BQ>
zz+}ClE~3CrkgA>pJfxNf^=-$BUD*H%(CW-OzAfD?EuWtUK_G=fs@}N|<_x9Hi!zS2
z11UlJ0X-R}qma?_bSF$C*K19ld2SlB&SNIVBvAal;tW_@36ZPgtc%rp+(ejw9_Fq)
zV6G98_o%L7WDbG~T#7U7WjKPa^F2q>E4#8iU_W`^+&72{wvZ-91>JugQt|N2(QdJP
z{-nv}rv>!6q}e^Vss}5hb^iVj9!Ri+($dj|y?hB!d++U6^>6=djxLfIRGS;>HHd1i
zZ&Qo&ok<8^a*iizXA4`_Cbr)&G&C$76wvC~@$qwU5lPandPxrzIeneR8OMgUOf^`Y
zXmO=i?yiMVEA3xY@UuuY#Z=D&dMBPhjN%@v3Zs>Eb@Fm@9JE2sqt!SmO)Kx;bI%_EZ<-BWpD8Q@z$IdxLx
zWF;j%4)%5c5S`C0l=Iv3aCKb+(6Jz$>}|_$>%J>_*;kDiCQ^ehkr40vg^c&XAo=C@
zR>_zd?uY*NUovVKj0;aaf75cOq+GGXA9HkzJFXG4HQ$PRsLe!a(H{P3knKulp4opZS3rgj|WsAhpK0$FR9Cz
zkgv?)1G2NTO9vPC_RK;*`uO|%`}hFD$Q;O=y#7lvI=zo18RiLE9vvOUVq-Sz3(PDn
zRrNlKd|8eD{eWti09v&OtfTlP6|-a_P0o>IUw!oZ%kM7hF?+KyuBSX
z0eNMn4YEZ-mS^&uLi617wNm@#Z
zuyS~?LP^Kr12{JLwSY8gRHWG->6V!>c^NG&xc;j?X2tcrk;4
zN5bll0dZnn+A{Z0-C17j@{@k9!buSdNM|P=={lt-9AQWjP?vK+_U(!+$j=6Yct-2V
z0Bt`57xcsnpgp=oh~=IuB!$92mf)p{!0BEx(g0N;L-Z5xWds-5dUQuYHn2PeCUy96
zIb2kWgCDap>FkNu?ouA#OB^KUt7cl4_iulTM8ocE7vaEhgyL
zg1j$6%bt8u^m8tV+Cnq9Cex=Ux#We5+FqP`a*sr4^|!DHD2FSXu%BN14f0U;GI9me
zUeAW3;?Ouf7)FzI1lMyf<`!3+42_e-{Z#XRmxXHl_A3{jL6H2ti)hi#6w)(5UucUH
z1{x1Y7deV+^tj2T01j1`4TftrZl=$jgv2fN+|H$a1{|KO=lb`(y}g|cV0Uwar8j41
zXTc8k-6z1oEW3Sey5V)#TdhZom;-E{as^rTr)To1iW9cBw!o%HmiI5vsjjVknRZJq
zrRNqq`)`{7GN}-Uexw>YPDQ2Zb-;*NN{4ECm!mjqs7hwfkF709;ZBC%kgs*(yLH&S
z9)0V>WIJXW40Mr=7W@=Y^U|9gn>b*Z>#>_=O1;d{j@I&)H+gB`KPddI13rZv+xn*%ArEYqwQ-I0>
zRODN0a7bblEaeIf@I=B=Q~0t-0|d{bNxeciPD54+V23B4FC-0EI)f%9Xciv%vQgb*
z1f|XrKi^$z;e32%V0btr?Uu9i%J=WKQc3sO>7mW=LIgvSKdBZVu5rSav=onXb8`nB
zW3zBRa`KXrXMjB}_-J8H+f~vYcliopSF()fIuJzsQ+&4u0kz*+Buace**ZXS7}@yDJ{L+Mg`I-~O!Bv}u{I4H
z?dC6Lwa^KrI)g`*1D*%mB`YUaTU!epU4MT+NG?vbgsG~kng>xpzEgIWWK9@J+qJTL
zDnLB&KlGDJp$h`>rqO>V&VPdXkvY(?ox{YdMQttTKNB%@
zT@+xta(AdYkKN(@`@>(r8t%#iC4`?|NK}QaPSRk3iYfrs<9;-~fjc`nv`L9zD60Ad
zcImlaNpmpbB$ew7(7=^Hd6=RFzC2Rsg3nSmW|Ko9Te^N|?(K0%{JVpG6PB~+*#{Uj
zzy|VA6u^BHUH;XGnGObILrU@%_yUp?Kx2*{@L6cZ6Hl(+`i?5gAvgKZ`V_zI@tvp$
zfFKW?fF8X&J86+A=qRB4eFBmj5jpGG4}=RyZ1R?VjNPEED59rqK_MVX3AC2i(&dNw
zxh>&X|8`*l@-zk6~d@hV~7Ds6VyX5m8Noz0n#0(#`l3d
z|AHj=4$~Cr|z8CM7E?J2E`X
zN=*TbkC&IX2UoU`3f12Rh)|Ijzr3ual4^Z+?usXv|J2?9Kl}-z6fDRgc(u>gLsjtr
z13=}!-Y5m?2Z&ssn12AI%g)IO>;rIqH+FWjKg8id719XU#x3m6h>HbL5gUExXzk&_
z?oUgXqM>>D$4yn4!=Sm&KbGFli)n#=O0s$nx<~@5e4MV4r*i+{6b3a|WnE~T1gvY;
zzv~CNg|$E*Oll^Fw856QlKL+zFSMv3Y_HjXpc&8u!qO`ME8SdOfu~F7V+mE!qnGvB
z26A@V!aFb{7IM4w=UMCfY>k^;W{f_!kTw3a>MA5(((=?>$ZtaKiM|{MlU%czFfkIP
z9v=VD_i5t01AZ(`6U3?gmJ@9XO2ayLBDo@}kDIT-(WFp)C7^G7k-c0;N%8uA8b=fX
zGo@r^5~OY{ExojJxX~ehxVO=fSs)X94si$6*M$V9sl{fU2Ij-oa{sYTeO_8U15R)ggC6gzHn
zRYiuWd7>h@T-n;%KzGLLwxT$D!{Y5Ij-Oqdko2I}hF@kAWTU_KZbg@r@|lX
zHUyT7N5R}Lq`#LUTn&r>$=5r%_ODoEJz1zJM5VNt0zrw((MH|o+{@C^*BX0sCT?ZS
zc4S^a9p)an5S-R(n|7MppD{$4UfF+diZfH=f$SArnmNHqqR)a&MAhjT7n|r!+!m%*
zvrzuX3aEU=VbO_95)Q(`H@lu`v_+QJy2rPML5^zC@_C
z@`L5xtp7oFF^dh-O2};6bvz$U$4NRmor?-lBKf#*XvhygZ-BN`=FpirUTfE
z27YjCfB8XVI3nUXOV*7{@6u%EPetPMnd-vdr!5cG;~K2a>?Cp|clO_*ew99he52ld
zUgHvJ0c`U)$gl0{3FC{sxU7uR_fWapgoVOTZjD1cf&RJi*3Z(Gy`2vK*ot0x6(|#p
z!d}Ils(TljfIOr*&9N@#m|7^wJu!7rm!%X&h3G*>28dAOaUk!Y{v25FS*nQV;QH{{
zUY%P<&c&D>v6wJ5Dej)o@-|0HpS31NOP`BO)I>g{q>P_v3=-4A$}|s6NebOji_a1@
z@5^}<7nTeEtlE2He)rA01=}UBdV5A13Nlm?FNyy_i2rVde-wP)AKPMgti%YV{sR;g
zq2IA;kmVjT5;ub
zPF$T$fdo>|>q4mnAOqHu{t8duv@q3kjGPdp@NOUr)+Smr>w&wF*YgY#P%VW)(qEe*
zVYnO%e0$974iDKHh(eAJL*y?!KBLxybWb>YIQJB_$3l_VgnOUhvo}hN%y%m`bY#;b
zNgluL6J`eecztW{$$Uew73JSjZ|b0`w@Z8jQh&lDM5p#Z14ECRHgwG45bIMO2Z_pA
zkML|&z3AEZ9K73*aJTxGF2$>tj*?<0vfYl@Zyq9E3y!Du^z?w#sh{7$
zcmvsG1Yl8NVPQdy6A2Q|02nR(2^T+6t%f&pRuf_?`#Xc?$v=z33IB(j~q+rjZy&2BKhwopUzH?`@>^zDY>(2^^#Zu;MsM{&D%RS
z>yvH(5%(HocahYKOC}*nRzHChVN+An{(g#A5+dz1MgW
z^QAwu2oDjL3Rp}2l;g-rIlnHVIeD&bazbORmjcXHNuXhcKltc(>EMtu&6xmp#h;%x
z@S0N=Wp#7b%CrEyRXORPuY7RCh{XZH4
zEndNM2GHt~9ZmHo=CUCHmkNS#o9GU2aN5t>N}~^58J}|)Q^`TOj7wq&)?Z(&d|Me?
z>IYS;fXyO1Wvw4G>&ezKpLrk;W+K5%6Z1htGcscVXljtbiFadW2-wk0NkfdNi$B>z
z^V3H-%$^oShmEBZ*~DGR_eT1ew#Eiuaz3UNzqAH5CRI%%8L;5AH4nCyi_cGjK2t&v
zF%}BSQ?S%A;ygLOGfHPWqR148+34tIW}{}U_mWKau5S_UBCgow#KeU?P|?5_{6Z8K
z3`Dkow#2rcD58?qglq2cFL+90btq_mR(NQ&#cpwx>Va|pSGjbfn9pOsU;sz!gVrVA-N0P-P*
zf&VC9+8)v&Qkwnq?HI`j5|7XQ8NqEr;K`2v&-H_XI?4GU|G&+lUfDDG0Xr#_7npH8
zK34iTfArO_>7pMUNh2n%vph|M5eh
zXF!`@{y<6zDT9HrvABv#+7nPLARFUPCR`+0Lhm>`hyLWY2CA+k3fH($keRE82aw5J
z{rwBG0DH2uu;Ad}xX;1}i2?Re&j1(_NSU{^w17^#<3Xm`@L&7e2vl1KT=*xZq^9_T
zDJR#aT8PxD?U$wt)x_BiddSRZ{&@OZUYgSfm{hVyc5*6i`d58ZiW_@#y!x)r?9VzG
z03gr#101jY2GkX+;>-g`19pGyj#VD?XTWPt01JWbL?D&*t`gVVY{m<)6l40IOLa7W
zwlnW|fGWs9=y)A2NF4`OIGrd2PJi4gHtO^S&AOxj$^(z@KN=ndsy75zL2oG1`%lDj
zC!;CH<+8zCLXgidPx6n_Xv)B}YmkI*K6euWI9q0A-;h=CO?Z*;}JJAI{SO+foSO$<&{rKoA!$gggN6jC$
z<64c`S`1=w@gmf_Ww)0(3mgEb8Osr&u6min;5#HXbzz|=3*F|}@wJ5KzP(MC_~k@1
z0knD~7qt>j+gJ4UM_0#8?V=yP!QcQ>=`!zk?|!;LKN%M>jxem?Wvw2Wj$H{yg^u~P
zX_E+huMIH2;;%Kof8cv>(*d%jxs`46v85;l?ibx?(Eq0CwL|Wagh-f;>GvOQ*(V3D
z)7Ne->h@~Ld6k;)@4lp;asUOF5{%)Q1XtOI=f4S}_45v`=3P_L^ytwmI&8=zvR(2r#oJ~5{&YTciTW>3OT#K9&OOaibzJ8OlFei>2
zS<{-6(WjHro$&bD-K(ewuJRj-CpK3oHJxik2TugAl9&>QmA}P*d?91u;lW7cuE}5F
zT+c^w0W$VhWuJdUaQ%?wkOq5)!f}umb%>rJ1v4g>3pBGBB=(hF!IXDI42(6LCi>>I
zSJW8e2JhG|ZqT
zf1@A^d|`bm3?118O-X*;sNO^<1YW0zW@pi5hI9*!*fNLSE;tR4M5h7wrHgOEy{OU;`i)SuYbCnRl^6R6D<)
zT@$;$SHkTBc5I8TyJhd8@IH+$zWcO$H?yS|_H4!bc(S~m%P3{-4uWriKW=zzmPV%{
zKB{*hh;U=RzXOwQVi$bBf@*+88WAI<1Hxu@)LO|V#3elq7L;=XS$uI(
za~-yC?(VN2f*BgOzOIQgZ5l@1zOJxA<<-6(qKD-|!Pg4oUFLph+w7>&+L}${&CSgo
zJ<
zyQ8m8b##HzH0I@STHMpQs>*i~eL`>e{n`j!xp&W8H>+GN{mc~;ySx#fVNbPN
zJ%PKCXL6_#+>&8(Wlc`D;y2N!lWmO)`Fx}I5Ihrs
z&KGW8UXoGMmU?H14TK!aUkd^ie3@OiL1nNKAYf(@nI+z9DL$RAt&Q>Vi8;KQp(I
ztyjR3p%lds2!CETppf1A6#Xu=DdhvK5dQh@3uT0;AzpU=BxWZUb0?(rAs9SaxlxXR
z*DI7=I0_Zp)N#k7kKS~OxyAG1PVTwP^9$I}Se^#0NX(?Cz1Z??w0gn8J8)LD+e^qs
zk06I8lKf&7JrTzYeGw4Tk-}8z)GapGE
zn?}>%Q$MrNxpE<_GmbXf2y38Z>udxr9eorZWEU^Lmo`?nA>+Kg8ER7dMgFxWSI5yt
z7a0y=^JgdP_OSk^XFW%4PrWepUemVZtd1n+9WI_@qoz?xDGaK
zs(W7kJR<-0XD1D87(S5UV?r6a+mNe5=_bIO@cHR;?TkjH@(w~&ObkZe>dP~b1fexnCi5FAv>
z!1^=A92l9Q%j&nE6~`Hx_%Q|QuZG*bt__dh4vF8}6HHR&(dLPo)Y+P339(aydB7ik
zOyNQ_cGH(D-~aGwLwH
zqvNLwqn%MrmU)&AC1jx};rthT`E@y!oYA9dm^<=>O@X%F8CF!V)>S>*Z^sNO72-I8
zn+glJtQahLuMNH>sNvEcstw4?S{cW4_Lsw}8f>l14>tNQR=-Zu-1R
zjT5#4U}^i2~=x(6ABUfc4cB1`)hXJcv*N*1m*1cyVtPqjM{a3
zI8i92pqR?Cf|V7_yPXS7mObrmz7)EZ7wpCJ2Qsf&mV>YM-kjLs}-
z7#Xvy0{7i8H9^U*H{B%D*d=%J+s6JFF_9Uye|g~;cKcIo6zqHkhOX}imC4Xr1I7t6
zPiM2wtGj<1HO+82o?uFoY8s6sxe)c0x9TP7a%3&Sd&iqXfuA1@%A-J0^fNBh
z$OEE+*hvP!>Hord|4_1i9
+)
+
+target_link_libraries(viewmode_check PMusr ${ROOT_LIBRARIES})
+
+add_test(
+  NAME viewmode-modes
+  COMMAND viewmode_check
+)
diff --git a/tests/viewmode_check/viewmode_check.cpp b/tests/viewmode_check/viewmode_check.cpp
new file mode 100644
index 00000000..10a59253
--- /dev/null
+++ b/tests/viewmode_check/viewmode_check.cpp
@@ -0,0 +1,88 @@
+/***************************************************************************
+
+  viewmode_check.cpp
+
+  Unit test of PMusrViewMode, the class which decouples the musrview plot
+  views (time/Fourier, difference d-f and f-d, average) from each other.
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+***************************************************************************/
+
+#include "PMusrViewMode.h"
+#include 
+#include 
+using namespace musrview;
+
+static void show(const char *what, const PMusrViewMode &m) {
+  std::cout << what << " -> " << m.GetKey() << "  |  " << m.GetLabel() << "\n";
+}
+
+int main() {
+  PMusrViewMode m;
+  show("initial", m);
+  assert(m.GetKey() == "time");
+
+  // d then f  => d-f
+  m.ToggleDifference(); show("d", m);
+  assert(m.GetKey() == "time+diff");
+  m.ToggleFourier();
+  show("d,f", m);
+  assert(m.IsDiffThenFourier());
+  assert(m.GetKey() == "fourier-power+diff-df");
+
+  // now 'd' again: old code refused (d-f-d), new one simply drops the difference
+  m.ToggleDifference(); show("d,f,d", m);
+  assert(m.GetKey() == "fourier-power");
+
+  // f then d => f-d
+  m.Reset();
+  m.ToggleFourier(); show("f", m);
+  m.ToggleDifference(); show("f,d", m);
+  assert(m.IsFourierThenDiff());
+  assert(m.GetKey() == "fourier-power+diff-fd");
+
+  // average is orthogonal: fourier-difference-average
+  m.ToggleAverage(); show("f,d,a", m);
+  assert(m.GetKey() == "fourier-power+diff-fd+avg");
+  assert(m.IsFourierThenDiff() && m.IsAverage() && m.IsFourier());
+
+  // component switch keeps difference + average
+  m.SetFourierComponent(EFourierComponent::kReal); show("+real", m);
+  assert(m.GetKey() == "fourier-real+diff-fd+avg");
+  assert(m.IsPhaseAdjustable());
+  m.SetFourierComponent(EFourierComponent::kPower);
+  assert(!m.IsPhaseAdjustable());
+
+  // back to time domain keeps difference + average
+  m.SetFourierComponent(EFourierComponent::kReal);
+  m.ToggleFourier(); show("back to time", m);
+  assert(m.GetKey() == "time+diff+avg");
+  assert(!m.IsPhaseAdjustable());
+
+  // and forward again -> difference was made in time domain => d-f
+  m.ToggleFourier(); show("forward", m);
+  assert(m.GetKey() == "fourier-real+diff-df+avg");
+
+  // capabilities
+  PMusrViewMode n;
+  n.SetFourierAllowed(false);
+  assert(!n.ToggleFourier());
+  assert(n.IsTimeDomain());
+  n.SetAverageAllowed(false);
+  assert(!n.ToggleAverage());
+  assert(!n.IsAverage());
+  show("non-musr", n);
+
+  // capability revoked while active
+  PMusrViewMode o;
+  o.ToggleFourier(); o.ToggleAverage();
+  o.SetFourierAllowed(false); o.SetAverageAllowed(false);
+  assert(o.IsTimeDomain() && !o.IsAverage());
+  show("revoked", o);
+
+  std::cout << "\nJSON: " << m.ToJson() << "\n";
+  std::cout << "\nall mode assertions passed\n";
+  return 0;
+}