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 00000000..6c7ea008
Binary files /dev/null and b/tests/musrview_check/ref/musrview-nonMusr-GapIntegrals/test-libGapIntegrals-ASCII_0.png differ
diff --git a/tests/viewmode_check/CMakeLists.txt b/tests/viewmode_check/CMakeLists.txt
new file mode 100644
index 00000000..ae563e70
--- /dev/null
+++ b/tests/viewmode_check/CMakeLists.txt
@@ -0,0 +1,15 @@
+#--- viewmode_check: unit test of the decoupled musrview view modes ----------
+
+add_executable(viewmode_check viewmode_check.cpp)
+
+target_include_directories(viewmode_check
+  BEFORE PRIVATE
+    $
+)
+
+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;
+}