From 2d2594c903ef619bd3dcd7499fac43c170ed58e0 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Fri, 13 Feb 2015 16:56:19 +0100 Subject: [PATCH] first almost feature complete version of musrFT --- ChangeLog | 8 +- src/classes/Makefile.am | 5 + src/classes/PFourierCanvas.cpp | 1312 +++++++++++++++++++++++++++ src/classes/PMusrCanvas.cpp | 203 +---- src/classes/PPrepFourier.cpp | 53 +- src/include/PFourier.h | 2 + src/include/PFourierCanvas.h | 161 ++++ src/include/PFourierCanvasLinkDef.h | 39 + src/include/PMusrCanvas.h | 4 +- src/include/PPrepFourier.h | 13 +- src/musrFT.cpp | 370 ++++++-- src/musrview.cpp | 5 +- 12 files changed, 1900 insertions(+), 275 deletions(-) create mode 100644 src/classes/PFourierCanvas.cpp create mode 100644 src/include/PFourierCanvas.h create mode 100644 src/include/PFourierCanvasLinkDef.h diff --git a/ChangeLog b/ChangeLog index 4f92c0c7..fcceb54f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,10 +4,10 @@ changes since 0.13.0 =================================== -NEW 2015-02-02 first partial implementation of a standalone Fourier transform/plotter: - musrFT. A LOT of the functionality is still missing! Initially it is - meant to be used for HAL-9500, i.e. Fourier transform WITHOUT lifetime - correction. +NEW 2015-02-13 first implementation of a standalone Fourier transform/plotter: + musrFT. Initially it is meant to be used for HAL-9500, + i.e. Fourier transform WITHOUT lifetime correction. + A first simple minded lifetime correction is implemented as well. NEW 2014-12-18 first implementation of a GLOBAL block which allows to shorten a typical msr-file. Duplicate entries from the RUN blocks can be added here. Furthermore, the 'lifetimecorrection' flag is diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am index 0be0160e..16ede5e0 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -4,6 +4,7 @@ h_sources = \ ../include/PFitterFcn.h \ ../include/PFitter.h \ ../include/PFourier.h \ + ../include/PFourierCanvas.h \ ../include/PFunctionGrammar.h \ ../include/PFunction.h \ ../include/PFunctionHandler.h \ @@ -27,6 +28,7 @@ h_sources_userFcn = \ ../include/PUserFcnBase.h h_linkdef = \ + ../include/PFourierCanvasLinkDef.h \ ../include/PMusrCanvasLinkDef.h \ ../include/PMusrT0LinkDef.h \ ../include/PStartupHandlerLinkDef.h @@ -35,6 +37,7 @@ h_linkdef_userFcn = \ ../include/PUserFcnBaseLinkDef.h dict_h_sources = \ + PFourierCanvasDict.h \ PMusrCanvasDict.h \ PMusrT0Dict.h \ PStartupHandlerDict.h @@ -46,6 +49,7 @@ cpp_sources = \ PFitter.cpp \ PFitterFcn.cpp \ PFourier.cpp \ + PFourierCanvas.cpp \ PFunction.cpp \ PFunctionHandler.cpp \ PMsr2Data.cpp \ @@ -68,6 +72,7 @@ cpp_sources_userFcn = \ PUserFcnBase.cpp dict_cpp_sources = \ + PFourierCanvasDict.cpp \ PMusrCanvasDict.cpp \ PMusrT0Dict.cpp \ PStartupHandlerDict.cpp diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp new file mode 100644 index 00000000..a8466779 --- /dev/null +++ b/src/classes/PFourierCanvas.cpp @@ -0,0 +1,1312 @@ +/*************************************************************************** + + PFourierCanvas.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2015 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 +#include +using namespace std; + +#include +#include +#include +#include +#include + +#include "PFourierCanvas.h" + +#define YINFO 0.2 +#define YTITLE 0.95 + +ClassImpQ(PFourierCanvas) + +//--------------------------------------------------------------------------- +/** + *

Constructor + */ +PFourierCanvas::PFourierCanvas() +{ + fTimeout = 0; + fTimeoutTimer = 0; + + fBatchMode = false; + fValid = false; + fAveragedView = false; + fCurrentPlotView = FOURIER_PLOT_NOT_GIVEN; + fInitialXRange[0] = 0.0; + fInitialXRange[1] = 0.0; + + fTitle = TString(""); + fXaxisTitle = TString(""); + + fCurrentFourierPhase = 0.0; + fCurrentFourierPhaseText = 0; + + fStyle = 0; + fImp = 0; + fBar = 0; + fPopupMain = 0; + fPopupSave = 0; + fPopupFourier = 0; + + fMainCanvas = 0; + fTitlePad = 0; + fFourierPad = 0; + fInfoPad = 0; +} + +//--------------------------------------------------------------------------- +/** + *

Constructor + * + * \param fourier + * \param title + * \param showAverage + * \param fourierPlotOpt + * \param fourierXrange + * \param phase + * \param wtopx + * \param wtopy + * \param ww + * \param wh + * \param batch + */ +PFourierCanvas::PFourierCanvas(vector &fourier, const Char_t* title, const Bool_t showAverage, + const Int_t fourierPlotOpt, Double_t fourierXrange[], Double_t phase, + Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, + const Bool_t batch) : + fBatchMode(batch), fAveragedView(showAverage), fCurrentPlotView(fourierPlotOpt), + fTitle(title), fFourier(fourier), fCurrentFourierPhase(phase) +{ + fInitialXRange[0] = fourierXrange[0]; + fInitialXRange[1] = fourierXrange[1]; + + fTimeout = 0; + fTimeoutTimer = 0; + + fValid = false; + + fCurrentFourierPhaseText = 0; + + CreateXaxisTitle(); + CreateStyle(); + InitFourierDataSets(); + InitFourierCanvas(title, wtopx, wtopy, ww, wh); + + gStyle->SetHistMinimumZero(kTRUE); // needed to enforce proper bar option handling +} + +//--------------------------------------------------------------------------- +/** + *

Constructor + * + * \param fourier + * \param title + * \param showAverage + * \param fourierPlotOpt + * \param fourierXrange + * \param phase + * \param wtopx + * \param wtopy + * \param ww + * \param wh + * \param markerList + * \param colorList + * \param batch + */ +PFourierCanvas::PFourierCanvas(vector &fourier, const Char_t* title, const Bool_t showAverage, + const Int_t fourierPlotOpt, Double_t fourierXrange[], Double_t phase, + Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, + const PIntVector markerList, const PIntVector colorList, const Bool_t batch) : + fBatchMode(batch), fAveragedView(showAverage), fCurrentPlotView(fourierPlotOpt), + fTitle(title), fFourier(fourier), fCurrentFourierPhase(phase), + fMarkerList(markerList), fColorList(colorList) +{ + fInitialXRange[0] = fourierXrange[0]; + fInitialXRange[1] = fourierXrange[1]; + + fTimeout = 0; + fTimeoutTimer = 0; + + fValid = false; + + fCurrentFourierPhaseText = 0; + + CreateXaxisTitle(); + CreateStyle(); + InitFourierDataSets(); + InitFourierCanvas(title, wtopx, wtopy, ww, wh); + + gStyle->SetHistMinimumZero(kTRUE); // needed to enforce proper bar option handling +} + +//--------------------------------------------------------------------------- +/** + *

Destructor + */ +PFourierCanvas::~PFourierCanvas() +{ + cout << endl << "debug> in ~PFourierCanvas()." << endl; + + if (fTimeoutTimer) + delete fTimeoutTimer; + + if (fCurrentFourierPhaseText) + delete fCurrentFourierPhaseText; + +/* + if (fStyle) { + delete fStyle; + fStyle = 0; + } +*/ + if (fTitlePad) { + fTitlePad->Clear(); + delete fTitlePad; + fTitlePad = 0; + } + + if (fFourierHistos.size() > 0) { + for (unsigned int i=0; iClear(); + delete fFourierPad; + fFourierPad = 0; + } +*/ + if (fInfoPad) { + fInfoPad->Clear(); + delete fInfoPad; + fInfoPad = 0; + } + +/* + if (fMainCanvas) { + delete fMainCanvas; + fMainCanvas = 0; + } +*/ +} + +//-------------------------------------------------------------------------- +// Done (SIGNAL) +//-------------------------------------------------------------------------- +/** + *

Signal emitted that the user wants to terminate the application. + * + * \param status Status info + */ +void PFourierCanvas::Done(Int_t status) +{ + Emit("Done(Int_t)", status); +} + +//-------------------------------------------------------------------------- +// HandleCmdKey (SLOT) +//-------------------------------------------------------------------------- +/** + *

Filters keyboard events, and if they are a command key (see below) carries out the + * necessary actions. + *

Currently implemented command keys: + * - 'q' quit musrview + * - 'u' unzoom to the original plot range given in the msr-file. + * - 'a' toggle between average view and single Fourier histo view. + * - '+' increment the phase (real/imaginary Fourier). The phase step is defined in the musrfit_startup.xml + * - '-' decrement the phase (real/imaginary Fourier). The phase step is defined in the musrfit_startup.xml + * + * \param event event type + * \param x character key + * \param y not used + * \param selected not used + */ +void PFourierCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected) +{ + if (event != kKeyPress) + return; + + if (fBatchMode) + return; + + if (x == 'q') { // quit + Done(0); + } else if (x == 'u') { // unzoom + if (fAveragedView) + PlotAverage(); + else + PlotFourier(); + } else if (x == 'a') { // toggle between average view and single histo view + // toggle average view flag + fAveragedView = !fAveragedView; + // update menu + if (fAveragedView) { + fPopupMain->CheckEntry(P_MENU_ID_AVERAGE); + HandleAverage(); + PlotAverage(); + } else { + fPopupMain->UnCheckEntry(P_MENU_ID_AVERAGE); + CleanupAverage(); + PlotFourier(); + } + } else if (x == '+') { // increment phase (Fourier real/imag) + IncrementFourierPhase(); + } else if (x == '-') { // decrement phase (Fourier real/imag) + DecrementFourierPhase(); + } +} + +//-------------------------------------------------------------------------- +// HandleMenuPopup (SLOT) +//-------------------------------------------------------------------------- +/** + *

Handles the MusrFT menu. + * + * \param id identification key of the selected menu + */ +void PFourierCanvas::HandleMenuPopup(Int_t id) +{ + if (fBatchMode) + return; + + if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL) { + // uncheck fourier popup items + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL); + fCurrentPlotView = FOURIER_PLOT_REAL; + PlotFourier(); + } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG) { + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG); + fCurrentPlotView = FOURIER_PLOT_IMAG; + PlotFourier(); + } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG) { + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG); + fCurrentPlotView = P_MENU_ID_FOURIER_REAL_AND_IMAG; + PlotFourier(); + } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR) { + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR); + fCurrentPlotView = FOURIER_PLOT_POWER; + PlotFourier(); + } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE) { + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE); + fCurrentPlotView = FOURIER_PLOT_PHASE; + PlotFourier(); + } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_PLUS) { + IncrementFourierPhase(); + } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_MINUS) { + DecrementFourierPhase(); + } else if (id == P_MENU_ID_AVERAGE) { + // toggle average view flag + fAveragedView = !fAveragedView; + // update menu + if (fAveragedView) { + fPopupMain->CheckEntry(P_MENU_ID_AVERAGE); + HandleAverage(); + PlotAverage(); + } else { + fPopupMain->UnCheckEntry(P_MENU_ID_AVERAGE); + CleanupAverage(); + PlotFourier(); + } + } else if (id == P_MENU_ID_SAVE_DATA+P_MENU_ID_SAVE_ASCII) { + cout << endl << "debug> Save still missing ..." << endl; + } +} + +//-------------------------------------------------------------------------- +// LastCanvasClosed (SLOT) +//-------------------------------------------------------------------------- +/** + *

Slot called when the last canvas has been closed. Will emit Done(0) which will + * terminate the application. + */ +void PFourierCanvas::LastCanvasClosed() +{ + if (gROOT->GetListOfCanvases()->IsEmpty()) { + Done(0); + } +} + +//-------------------------------------------------------------------------- +// UpdateFourierPad (public) +//-------------------------------------------------------------------------- +/** + *

Feeds the pad with the Fourier data sets and refreshes it. + */ +void PFourierCanvas::UpdateFourierPad() +{ + if (!fValid) + return; + + if (fAveragedView) + PlotAverage(); + else + PlotFourier(); + + fFourierPad->Draw(); + fMainCanvas->cd(); + fMainCanvas->Update(); +} + +//-------------------------------------------------------------------------- +// UpdateInfoPad (public) +//-------------------------------------------------------------------------- +/** + *

Feeds the pad with the legend and refreshes it. + */ +void PFourierCanvas::UpdateInfoPad() +{ + if (!fValid) + return; + + // write header line + TDatime dt; + char dtStr[128]; + strncpy(dtStr, dt.AsSQLString(), sizeof(dtStr)); + TString str("musrFT: "); + str += dtStr; + fInfoPad->SetHeader(str); + + char title[1024]; + for (unsigned int i=0; iGetDataTitle(), sizeof(title)); + // add entry + fInfoPad->AddEntry(fFourierHistos[i].dataFourierRe, title, "p"); + } + + fInfoPad->Draw(); + fMainCanvas->cd(); + fMainCanvas->Update(); +} + +//-------------------------------------------------------------------------- +// SetTimeout (public) +//-------------------------------------------------------------------------- +/** + *

sets the timeout after which the program shall terminate. + * + * \param timeout after which the done signal shall be emitted. Given in seconds + */ +void PFourierCanvas::SetTimeout(Int_t timeout) +{ + fTimeout = timeout; + + if (fTimeout <= 0) + return; + + if (fTimeoutTimer) { + delete fTimeoutTimer; + fTimeoutTimer = 0; + } + fTimeoutTimer = new TTimer(); + + fTimeoutTimer->Connect("Timeout()", "PFourierCanvas", this, "Done()"); + + fTimeoutTimer->Start(1000*fTimeout, kTRUE); +} + +//-------------------------------------------------------------------------- +// SaveGraphicsAndQuit +//-------------------------------------------------------------------------- +/** + *

Will save the canvas as graphics output. Needed in the batch mode of musrview. + * + * \param fileName file name under which the canvas shall be saved. + */ +void PFourierCanvas::SaveGraphicsAndQuit(const Char_t *fileName) +{ + cout << endl << ">> SaveGraphicsAndQuit: will dump the canvas into a graphics output file: " << fileName << " ..." << endl; + + fMainCanvas->SaveAs(fileName); + + Done(0); +} + +//-------------------------------------------------------------------------- +// SaveDataAscii +//-------------------------------------------------------------------------- +/** + *

Saves the currently displayed Fourier data set in ascii column format. + */ +void PFourierCanvas::SaveDataAscii() +{ + // STILL MISSING +} + +//-------------------------------------------------------------------------- +// CreateXaxisTitle (private) +//-------------------------------------------------------------------------- +/** + *

Creates the x-axis title based on the Fourier units used. + */ +void PFourierCanvas::CreateXaxisTitle() +{ + switch (fFourier[0]->GetUnitTag()) { + case FOURIER_UNIT_GAUSS: + fXaxisTitle = TString("Field (Gauss)"); + break; + case FOURIER_UNIT_TESLA: + fXaxisTitle = TString("Field (Tesla)"); + break; + case FOURIER_UNIT_FREQ: + fXaxisTitle = TString("Frequency (MHz)"); + break; + case FOURIER_UNIT_CYCLES: + fXaxisTitle = TString("Angular Frequency (Mc/s)"); + break; + default: + fXaxisTitle = TString("??"); + break; + } +} + +//-------------------------------------------------------------------------- +// CreateStyle (private) +//-------------------------------------------------------------------------- +/** + *

Set styles for the canvas. Perhaps one could transfer them to the startup-file in the future. + */ +void PFourierCanvas::CreateStyle() +{ + TString musrFTStyle("musrFTStyle"); + fStyle = new TStyle(musrFTStyle, musrFTStyle); + fStyle->SetOptStat(0); // no statistics options + fStyle->SetOptTitle(0); // no title + fStyle->cd(); +} + +//-------------------------------------------------------------------------- +// InitFourierDataSets (private) +//-------------------------------------------------------------------------- +/** + *

Initialize the Fourier data sets, i.e. get the TH1F objects, set markers, + * set colors, ... + */ +void PFourierCanvas::InitFourierDataSets() +{ + // get all the Fourier histos + fFourierHistos.resize(fFourier.size()); + for (unsigned int i=0; iGetRealFourier(); + fFourierHistos[i].dataFourierIm = fFourier[i]->GetImaginaryFourier(); + fFourierHistos[i].dataFourierPwr = fFourier[i]->GetPowerFourier(); + fFourierHistos[i].dataFourierPhase = fFourier[i]->GetPhaseFourier(); + } + + // rescale histo to abs(maximum) == 1 + Double_t max = 0.0, dval = 0.0; + Int_t start = fFourierHistos[0].dataFourierRe->FindBin(fInitialXRange[0]); + Int_t end = fFourierHistos[0].dataFourierRe->FindBin(fInitialXRange[1]); + // real + for (unsigned int i=0; iGetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; iGetNbinsX(); j++) { + fFourierHistos[i].dataFourierRe->SetBinContent(j, fFourierHistos[i].dataFourierRe->GetBinContent(j)/fabs(max)); + } + } + // imaginary + for (unsigned int i=0; iGetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; iGetNbinsX(); j++) { + fFourierHistos[i].dataFourierIm->SetBinContent(j, fFourierHistos[i].dataFourierIm->GetBinContent(j)/fabs(max)); + } + } + // power + max = 0.0; + for (unsigned int i=0; iGetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; iGetNbinsX(); j++) { + fFourierHistos[i].dataFourierPwr->SetBinContent(j, fFourierHistos[i].dataFourierPwr->GetBinContent(j)/fabs(max)); + } + } + // phase + max = 0.0; + for (unsigned int i=0; iGetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; iGetNbinsX(); j++) { + fFourierHistos[i].dataFourierPhase->SetBinContent(j, fFourierHistos[i].dataFourierPhase->GetBinContent(j)/fabs(max)); + } + } + + // set the marker and line color + for (unsigned int i=0; iSetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierRe->SetLineColor(fColorList[i]); + fFourierHistos[i].dataFourierIm->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierIm->SetLineColor(fColorList[i]); + fFourierHistos[i].dataFourierPwr->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierPwr->SetLineColor(fColorList[i]); + fFourierHistos[i].dataFourierPhase->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierPhase->SetLineColor(fColorList[i]); + } else { + TRandom rand(i); + Int_t color = TColor::GetColor((Int_t)rand.Integer(255), (Int_t)rand.Integer(255), (Int_t)rand.Integer(255)); + fFourierHistos[i].dataFourierRe->SetMarkerColor(color); + fFourierHistos[i].dataFourierRe->SetLineColor(color); + fFourierHistos[i].dataFourierIm->SetMarkerColor(color); + fFourierHistos[i].dataFourierIm->SetLineColor(color); + fFourierHistos[i].dataFourierPwr->SetMarkerColor(color); + fFourierHistos[i].dataFourierPwr->SetLineColor(color); + fFourierHistos[i].dataFourierPhase->SetMarkerColor(color); + fFourierHistos[i].dataFourierPhase->SetLineColor(color); + } + } + + // set the marker symbol and size + Int_t style; + for (unsigned int i=0; iSetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierRe->SetMarkerSize(0.7); + fFourierHistos[i].dataFourierIm->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierIm->SetMarkerSize(0.7); + fFourierHistos[i].dataFourierPwr->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierPwr->SetMarkerSize(0.7); + fFourierHistos[i].dataFourierPhase->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierPhase->SetMarkerSize(0.7); + } else { + TRandom rand(i); + style = 20+(Int_t)rand.Integer(10); + fFourierHistos[i].dataFourierRe->SetMarkerStyle(style); + fFourierHistos[i].dataFourierRe->SetMarkerSize(0.7); + fFourierHistos[i].dataFourierIm->SetMarkerStyle(style); + fFourierHistos[i].dataFourierIm->SetMarkerSize(0.7); + fFourierHistos[i].dataFourierPwr->SetMarkerStyle(style); + fFourierHistos[i].dataFourierPwr->SetMarkerSize(0.7); + fFourierHistos[i].dataFourierPhase->SetMarkerStyle(style); + fFourierHistos[i].dataFourierPhase->SetMarkerSize(0.7); + } + } + + // initialize average histos + fFourierAverage.dataFourierRe = 0; + fFourierAverage.dataFourierIm = 0; + fFourierAverage.dataFourierPwr = 0; + fFourierAverage.dataFourierPhase = 0; +} + +//-------------------------------------------------------------------------- +// InitFourierCanvas (private) +//-------------------------------------------------------------------------- +/** + *

Initialize the class, and sets up the necessary objects. + * + * \param title Title to be displayed + * \param wtopx top x coordinate (in pixels) to place the canvas. + * \param wtopy top y coordinate (in pixels) to place the canvas. + * \param ww width (in pixels) of the canvas. + * \param wh height (in pixels) of the canvas. + */ +void PFourierCanvas::InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh) +{ + fImp = 0; + fBar = 0; + fPopupMain = 0; + fPopupSave = 0; + fPopupFourier = 0; + + fMainCanvas = 0; + fTitlePad = 0; + fFourierPad = 0; + fInfoPad = 0; + + // invoke canvas + TString canvasName = TString("fMainCanvas"); + fMainCanvas = new TCanvas(canvasName.Data(), title, wtopx, wtopy, ww, wh); + if (fMainCanvas == 0) { + cerr << endl << "PFourierCanvas::PFourierCanvas: **PANIC ERROR**: Couldn't invoke " << canvasName.Data(); + cerr << endl; + return; + } + + // add canvas menu if not in batch mode + if (!fBatchMode) { + fImp = (TRootCanvas*)fMainCanvas->GetCanvasImp(); + fBar = fImp->GetMenuBar(); + fPopupMain = fBar->AddPopup("&MusrFT"); + + fPopupFourier = new TGPopupMenu(); + + fPopupMain->AddPopup("&Fourier", fPopupFourier); + fPopupFourier->AddEntry("Show Real", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL); + fPopupFourier->AddEntry("Show Imag", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG); + fPopupFourier->AddEntry("Show Real+Imag", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG); + fPopupFourier->AddEntry("Show Power", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR); + fPopupFourier->AddEntry("Show Phase", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE); + fPopupFourier->AddSeparator(); + fPopupFourier->AddEntry("Phase +", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_PLUS); + fPopupFourier->AddEntry("Phase -", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_MINUS); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_PLUS); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_MINUS); + + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL); + break; + case FOURIER_PLOT_IMAG: + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG); + break; + case FOURIER_PLOT_REAL_AND_IMAG: + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG); + break; + case FOURIER_PLOT_POWER: + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR); + break; + case FOURIER_PLOT_PHASE: + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE); + break; + default: + break; + } + + fPopupMain->AddEntry("Average", P_MENU_ID_AVERAGE); + if (fAveragedView) + fPopupMain->CheckEntry(P_MENU_ID_AVERAGE); + fPopupMain->AddSeparator(); + + fPopupSave = new TGPopupMenu(); + fPopupSave->AddEntry("Save ascii", P_MENU_ID_SAVE_DATA+P_MENU_ID_SAVE_ASCII); + + fPopupMain->AddPopup("&Save Data", fPopupSave); + fBar->MapSubwindows(); + fBar->Layout(); + + fPopupMain->Connect("TGPopupMenu", "Activated(Int_t)", "PFourierCanvas", this, "HandleMenuPopup(Int_t)"); + } + + // divide the canvas into sub pads + // title pad + fTitlePad = new TPaveText(0.0, YTITLE, 1.0, 1.0, "NDC"); + if (fTitlePad == 0) { + cerr << endl << "PFourierCanvas::PFourierCanvas: **PANIC ERROR**: Couldn't invoke fTitlePad"; + cerr << endl; + return; + } + fTitlePad->SetFillColor(TColor::GetColor(255,255,255)); + fTitlePad->SetTextAlign(12); // middle, left + fTitlePad->AddText(title); + fTitlePad->Draw(); + + // fourier pad + fFourierPad = new TPad("fourierPad", "fourierPad", 0.0, YINFO, 1.0, YTITLE); + if (fFourierPad == 0) { + cerr << endl << "PFourierCanvas::PFourierCanvas: **PANIC ERROR**: Couldn't invoke fFourierPad"; + cerr << endl; + return; + } + fFourierPad->SetFillColor(TColor::GetColor(255,255,255)); + fFourierPad->Draw(); + + // info pad + fInfoPad = new TLegend(0.0, 0.0, 1.0, YINFO, "NDC"); + if (fInfoPad == 0) { + cerr << endl << "PFourierCanvas::PFourierCanvas: **PANIC ERROR**: Couldn't invoke fInfoPad"; + cerr << endl; + return; + } + fInfoPad->SetFillColor(TColor::GetColor(255,255,255)); + fInfoPad->SetTextAlign(12); // middle, left + + fValid = true; + + fMainCanvas->cd(); + + fMainCanvas->Show(); + + fMainCanvas->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "PFourierCanvas", + this, "HandleCmdKey(Int_t,Int_t,Int_t,TObject*)"); +} + +//-------------------------------------------------------------------------- +// CleanupAverage (private) +//-------------------------------------------------------------------------- +/** + *

Cleanup average Fourier data sets. + */ +void PFourierCanvas::CleanupAverage() +{ + if (fFourierAverage.dataFourierRe != 0) { + delete fFourierAverage.dataFourierRe; + fFourierAverage.dataFourierRe = 0; + } + if (fFourierAverage.dataFourierIm) { + delete fFourierAverage.dataFourierIm; + fFourierAverage.dataFourierIm = 0; + } + if (fFourierAverage.dataFourierPwr) { + delete fFourierAverage.dataFourierPwr; + fFourierAverage.dataFourierPwr = 0; + } + if (fFourierAverage.dataFourierPhase) { + delete fFourierAverage.dataFourierPhase; + fFourierAverage.dataFourierPhase = 0; + } +} + +//-------------------------------------------------------------------------- +// HandleAverage (private) +//-------------------------------------------------------------------------- +/** + *

Average Fourier data sets. + */ +void PFourierCanvas::HandleAverage() +{ + if (fFourierHistos.size() == 0) + return; + + CleanupAverage(); + + TString name(""); + + // create average histograms + name = TString(fFourierHistos[0].dataFourierRe->GetTitle()) + "_avg"; + fFourierAverage.dataFourierRe = new TH1F(name, name, fFourierHistos[0].dataFourierRe->GetNbinsX(), + fFourierHistos[0].dataFourierRe->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierRe->GetXaxis()->GetXmax()); + + name = TString(fFourierHistos[0].dataFourierIm->GetTitle()) + "_avg"; + fFourierAverage.dataFourierIm = new TH1F(name, name, fFourierHistos[0].dataFourierIm->GetNbinsX(), + fFourierHistos[0].dataFourierIm->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierIm->GetXaxis()->GetXmax()); + + name = TString(fFourierHistos[0].dataFourierPwr->GetTitle()) + "_avg"; + fFourierAverage.dataFourierPwr = new TH1F(name, name, fFourierHistos[0].dataFourierPwr->GetNbinsX(), + fFourierHistos[0].dataFourierPwr->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierPwr->GetXaxis()->GetXmax()); + + name = TString(fFourierHistos[0].dataFourierPhase->GetTitle()) + "_avg"; + fFourierAverage.dataFourierPhase = new TH1F(name, name, fFourierHistos[0].dataFourierPhase->GetNbinsX(), + fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmax()); + + Double_t dval=0.0; + // real average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; iGetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[i].dataFourierRe, fFourierHistos[0].dataFourierRe->GetBinCenter(j)); + } + fFourierAverage.dataFourierRe->SetBinContent(j, dval/fFourierHistos.size()); + } + // set marker color, line color, maker size, marker type + fFourierAverage.dataFourierRe->SetMarkerColor(kBlack); + fFourierAverage.dataFourierRe->SetLineColor(kBlack); + fFourierAverage.dataFourierRe->SetMarkerSize(0.8); + fFourierAverage.dataFourierRe->SetMarkerStyle(22); // closed up triangle + + // imaginary average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; iGetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[i].dataFourierIm, fFourierHistos[0].dataFourierIm->GetBinCenter(j)); + } + fFourierAverage.dataFourierIm->SetBinContent(j, dval/fFourierHistos.size()); + } + // set marker color, line color, maker size, marker type + fFourierAverage.dataFourierIm->SetMarkerColor(kBlack); + fFourierAverage.dataFourierIm->SetLineColor(kBlack); + fFourierAverage.dataFourierIm->SetMarkerSize(0.8); + fFourierAverage.dataFourierIm->SetMarkerStyle(22); // closed up triangle + + // power average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; iGetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[i].dataFourierPwr, fFourierHistos[0].dataFourierPwr->GetBinCenter(j)); + } + fFourierAverage.dataFourierPwr->SetBinContent(j, dval/fFourierHistos.size()); + } + // set marker color, line color, maker size, marker type + fFourierAverage.dataFourierPwr->SetMarkerColor(kBlack); + fFourierAverage.dataFourierPwr->SetLineColor(kBlack); + fFourierAverage.dataFourierPwr->SetMarkerSize(0.8); + fFourierAverage.dataFourierPwr->SetMarkerStyle(22); // closed up triangle + + // phase average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; iGetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[i].dataFourierPhase, fFourierHistos[0].dataFourierPhase->GetBinCenter(j)); + } + fFourierAverage.dataFourierPhase->SetBinContent(j, dval/fFourierHistos.size()); + } + // set marker color, line color, maker size, marker type + fFourierAverage.dataFourierPhase->SetMarkerColor(kBlack); + fFourierAverage.dataFourierPhase->SetLineColor(kBlack); + fFourierAverage.dataFourierPhase->SetMarkerSize(0.8); + fFourierAverage.dataFourierPhase->SetMarkerStyle(22); +} + +//-------------------------------------------------------------------------- +// PlotFourier (private) +//-------------------------------------------------------------------------- +/** + *

Plot the Fourier spectra. + */ +void PFourierCanvas::PlotFourier() +{ + fFourierPad->cd(); + + Double_t xmin=0, xmax=0; + xmin = fInitialXRange[0]; + xmax = fInitialXRange[1]; + + Double_t ymin=0, ymax=0; + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + fFourierHistos[0].dataFourierRe->GetXaxis()->SetRangeUser(xmin, xmax); + ymin = GetMinimum(fFourierHistos[0].dataFourierRe, xmin, xmax); + ymax = GetMaximum(fFourierHistos[0].dataFourierRe, xmin, xmax); + for (unsigned int i=1; i ymax) + ymax = GetMaximum(fFourierHistos[i].dataFourierRe, xmin, xmax); + if (GetMinimum(fFourierHistos[i].dataFourierRe, xmin, xmax) < ymin) + ymin = GetMinimum(fFourierHistos[i].dataFourierRe, xmin, xmax); + } + fFourierHistos[0].dataFourierRe->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + fFourierHistos[0].dataFourierRe->GetYaxis()->SetTitle("Real"); + fFourierHistos[0].dataFourierRe->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierHistos[0].dataFourierRe->Draw("p"); + for (unsigned int i=1; iDraw("psame"); + } + break; + case FOURIER_PLOT_IMAG: + fFourierHistos[0].dataFourierIm->GetXaxis()->SetRangeUser(xmin, xmax); + ymin = GetMinimum(fFourierHistos[0].dataFourierIm, xmin, xmax); + ymax = GetMaximum(fFourierHistos[0].dataFourierIm, xmin, xmax); + for (unsigned int i=1; i ymax) + ymax = GetMaximum(fFourierHistos[i].dataFourierIm, xmin, xmax); + if (GetMinimum(fFourierHistos[i].dataFourierIm, xmin, xmax) < ymin) + ymin = GetMinimum(fFourierHistos[i].dataFourierIm, xmin, xmax); + } + fFourierHistos[0].dataFourierIm->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + fFourierHistos[0].dataFourierIm->GetYaxis()->SetTitle("Imag"); + fFourierHistos[0].dataFourierIm->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierHistos[0].dataFourierIm->Draw("p"); + for (unsigned int i=1; iDraw("psame"); + } + break; + case FOURIER_PLOT_REAL_AND_IMAG: + fFourierHistos[0].dataFourierRe->GetXaxis()->SetRangeUser(xmin, xmax); + ymin = GetMinimum(fFourierHistos[0].dataFourierRe, xmin, xmax); + ymax = GetMaximum(fFourierHistos[0].dataFourierRe, xmin, xmax); + for (unsigned int i=1; i ymax) + ymax = GetMaximum(fFourierHistos[i].dataFourierRe, xmin, xmax); + if (GetMaximum(fFourierHistos[i].dataFourierIm, xmin, xmax) > ymax) + ymax = GetMaximum(fFourierHistos[i].dataFourierIm, xmin, xmax); + if (GetMinimum(fFourierHistos[i].dataFourierRe, xmin, xmax) < ymin) + ymin = GetMinimum(fFourierHistos[i].dataFourierRe, xmin, xmax); + if (GetMinimum(fFourierHistos[i].dataFourierIm, xmin, xmax) < ymin) + ymin = GetMinimum(fFourierHistos[i].dataFourierIm, xmin, xmax); + } + fFourierHistos[0].dataFourierRe->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + fFourierHistos[0].dataFourierRe->GetYaxis()->SetTitle("Real"); + fFourierHistos[0].dataFourierRe->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierHistos[0].dataFourierRe->Draw("p"); + for (unsigned int i=1; iDraw("psame"); + } + for (unsigned int i=0; iDraw("psame"); + } + break; + case FOURIER_PLOT_POWER: + // get maximum of Fourier data in the range + fFourierHistos[0].dataFourierPwr->GetXaxis()->SetRangeUser(xmin, xmax); + ymin = 0.0; + ymax = GetMaximum(fFourierHistos[0].dataFourierPwr, xmin, xmax); + for (unsigned int i=1; i ymax) + ymax = GetMaximum(fFourierHistos[i].dataFourierPwr, xmin, xmax); + } + fFourierHistos[0].dataFourierPwr->GetYaxis()->SetRangeUser(ymin, 1.03*ymax); + fFourierHistos[0].dataFourierPwr->GetYaxis()->SetTitle("Power"); + fFourierHistos[0].dataFourierPwr->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierHistos[0].dataFourierPwr->Draw("p"); + for (unsigned int i=1; iDraw("psame"); + } + break; + case FOURIER_PLOT_PHASE: + fFourierHistos[0].dataFourierPhase->GetXaxis()->SetRangeUser(xmin, xmax); + ymin = GetMinimum(fFourierHistos[0].dataFourierPhase, xmin, xmax); + ymax = GetMaximum(fFourierHistos[0].dataFourierPhase, xmin, xmax); + for (unsigned int i=1; i ymax) + ymax = GetMaximum(fFourierHistos[i].dataFourierPhase, xmin, xmax); + if (GetMinimum(fFourierHistos[i].dataFourierPhase, xmin, xmax) < ymin) + ymin = GetMinimum(fFourierHistos[i].dataFourierPhase, xmin, xmax); + } + fFourierHistos[0].dataFourierRe->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + fFourierHistos[0].dataFourierPhase->GetYaxis()->SetTitle("Phase"); + fFourierHistos[0].dataFourierPhase->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierHistos[0].dataFourierPhase->Draw("p"); + for (unsigned int i=1; iDraw("psame"); + } + break; + default: + break; + } + + fFourierPad->Update(); + + fMainCanvas->cd(); + fMainCanvas->Update(); +} + +//-------------------------------------------------------------------------- +// PlotFourierPhaseValue (private) +//-------------------------------------------------------------------------- +/** + *

Writes the Fourier phase value into the data window. + */ +void PFourierCanvas::PlotFourierPhaseValue() +{ + // check if phase TLatex object is present + if (fCurrentFourierPhaseText) { + delete fCurrentFourierPhaseText; + fCurrentFourierPhaseText = 0; + } + + double x, y; + TString str; + + // plot Fourier phase + str = TString("phase = "); + str += fCurrentFourierPhase; + x = 0.7; + y = 0.85; + fCurrentFourierPhaseText = new TLatex(); + fCurrentFourierPhaseText->SetNDC(kTRUE); + fCurrentFourierPhaseText->SetText(x, y, str.Data()); + fCurrentFourierPhaseText->SetTextFont(62); + fCurrentFourierPhaseText->SetTextSize(0.03); + + fFourierPad->cd(); + + fCurrentFourierPhaseText->Draw(); + + fFourierPad->Update(); +} + +//-------------------------------------------------------------------------- +// PlotAverage (private) +//-------------------------------------------------------------------------- +/** + *

Plot the average of the given Fourier data sets. + */ +void PFourierCanvas::PlotAverage() +{ + fFourierPad->cd(); + + if (fFourierAverage.dataFourierRe == 0) { // ups, HandlerAverage never called! + HandleAverage(); + } + + Double_t xmin=0, xmax=0; + xmin = fInitialXRange[0]; + xmax = fInitialXRange[1]; + + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + fFourierAverage.dataFourierRe->GetYaxis()->SetTitle(""); + fFourierAverage.dataFourierRe->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierAverage.dataFourierRe->GetXaxis()->SetRangeUser(xmin, xmax); + fFourierAverage.dataFourierRe->Draw("p"); + break; + case FOURIER_PLOT_IMAG: + fFourierAverage.dataFourierIm->GetYaxis()->SetTitle(""); + fFourierAverage.dataFourierIm->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierAverage.dataFourierIm->GetXaxis()->SetRangeUser(xmin, xmax); + fFourierAverage.dataFourierIm->Draw("p"); + break; + case FOURIER_PLOT_REAL_AND_IMAG: + fFourierAverage.dataFourierRe->GetYaxis()->SetTitle(""); + fFourierAverage.dataFourierRe->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierAverage.dataFourierRe->GetXaxis()->SetRangeUser(xmin, xmax); + fFourierAverage.dataFourierRe->Draw("p"); + fFourierAverage.dataFourierIm->Draw("psame"); + break; + case FOURIER_PLOT_POWER: + // get maximum of Fourier data in the range + fFourierAverage.dataFourierPwr->GetXaxis()->SetRangeUser(xmin, xmax); + fFourierAverage.dataFourierPwr->GetYaxis()->SetTitle(""); + fFourierAverage.dataFourierPwr->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierAverage.dataFourierPwr->Draw("p"); + break; + case FOURIER_PLOT_PHASE: + fFourierAverage.dataFourierPhase->GetYaxis()->SetTitle(""); + fFourierAverage.dataFourierPhase->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierAverage.dataFourierPhase->GetXaxis()->SetRangeUser(xmin, xmax); + fFourierAverage.dataFourierPhase->Draw("p"); + break; + default: + break; + } + + fFourierPad->Update(); + + fMainCanvas->cd(); + fMainCanvas->Update(); +} + +//-------------------------------------------------------------------------- +// IncrementFourierPhase (private) +//-------------------------------------------------------------------------- +/** + *

Increments the Fourier phase and recalculate the real/imaginary part of the Fourier transform. + */ +void PFourierCanvas::IncrementFourierPhase() +{ + if ((fCurrentPlotView == FOURIER_PLOT_POWER) || (fCurrentPlotView == FOURIER_PLOT_PHASE)) + return; + + double re, im; + double inc = 1.0; + const double cp = TMath::Cos(inc/180.0*TMath::Pi()); + const double sp = TMath::Sin(inc/180.0*TMath::Pi()); + + fCurrentFourierPhase += inc; + PlotFourierPhaseValue(); + + for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set + // calculate new fourier data set value + re = fFourierHistos[i].dataFourierRe->GetBinContent(j) * cp + fFourierHistos[i].dataFourierIm->GetBinContent(j) * sp; + im = fFourierHistos[i].dataFourierIm->GetBinContent(j) * cp - fFourierHistos[i].dataFourierRe->GetBinContent(j) * sp; + // overwrite fourier data set value + fFourierHistos[i].dataFourierRe->SetBinContent(j, re); + fFourierHistos[i].dataFourierIm->SetBinContent(j, im); + } + } + } +} + +//-------------------------------------------------------------------------- +// DecrementFourierPhase (private) +//-------------------------------------------------------------------------- +/** + *

Decrements the Fourier phase and recalculate the real/imaginary part of the Fourier transform. + */ +void PFourierCanvas::DecrementFourierPhase() +{ + if ((fCurrentPlotView == FOURIER_PLOT_POWER) || (fCurrentPlotView == FOURIER_PLOT_PHASE)) + return; + + double re, im; + double inc = 1.0; + const double cp = TMath::Cos(inc/180.0*TMath::Pi()); + const double sp = TMath::Sin(inc/180.0*TMath::Pi()); + + fCurrentFourierPhase -= inc; + PlotFourierPhaseValue(); + + for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set + // calculate new fourier data set value + re = fFourierHistos[i].dataFourierRe->GetBinContent(j) * cp + fFourierHistos[i].dataFourierIm->GetBinContent(j) * sp; + im = fFourierHistos[i].dataFourierIm->GetBinContent(j) * cp - fFourierHistos[i].dataFourierRe->GetBinContent(j) * sp; + // overwrite fourier data set value + fFourierHistos[i].dataFourierRe->SetBinContent(j, re); + fFourierHistos[i].dataFourierIm->SetBinContent(j, im); + } + } + } +} + +//-------------------------------------------------------------------------- +// 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 PFourierCanvas::GetMaximum(TH1F* histo, Double_t xmin, Double_t xmax) +{ + if (histo == 0) + 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 PFourierCanvas::GetMinimum(TH1F* histo, Double_t xmin, Double_t xmax) +{ + if (histo == 0) + 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; +} + +//-------------------------------------------------------------------------- +// 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 PFourierCanvas::GetInterpolatedValue(TH1F* histo, Double_t xVal) +{ + if (histo == 0) + 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; +} diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 9762e17b..4fcf2bd5 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -3797,15 +3797,10 @@ void PMusrCanvas::HandleAverage() // calculate all the average data sets double dval; if (fDataAvg.data != 0) { - if (!CalcAlignment(eTime)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: data: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].data->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].data, fData[0].data->GetBinContent(i)); } fDataAvg.data->SetBinContent(i, dval/fData.size()); } @@ -3816,15 +3811,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.data->SetMarkerStyle(fData[0].data->GetMarkerStyle()); } if (fDataAvg.dataFourierRe != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Re: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].dataFourierRe->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierRe, fData[0].dataFourierRe->GetBinContent(i)); } fDataAvg.dataFourierRe->SetBinContent(i, dval/fData.size()); } @@ -3835,15 +3825,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierRe->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle()); } if (fDataAvg.dataFourierIm != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Im: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].dataFourierIm->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierIm, fData[0].dataFourierIm->GetBinContent(i)); } fDataAvg.dataFourierIm->SetBinContent(i, dval/fData.size()); } @@ -3854,15 +3839,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierIm->SetMarkerStyle(fData[0].dataFourierIm->GetMarkerStyle()); } if (fDataAvg.dataFourierPwr != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Pwr: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].dataFourierPwr->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierPwr, fData[0].dataFourierPwr->GetBinContent(i)); } fDataAvg.dataFourierPwr->SetBinContent(i, dval/fData.size()); } @@ -3873,15 +3853,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierPwr->SetMarkerStyle(fData[0].dataFourierPwr->GetMarkerStyle()); } if (fDataAvg.dataFourierPhase != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Phase: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].dataFourierPhase->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierPhase, fData[0].dataFourierPhase->GetBinContent(i)); } fDataAvg.dataFourierPhase->SetBinContent(i, dval/fData.size()); } @@ -3892,15 +3867,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); } if (fDataAvg.theory != 0) { - if (!CalcAlignment(eTheoTime)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: theory: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].theory->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].theory, fData[0].theory->GetBinContent(i)); } fDataAvg.theory->SetBinContent(i, dval/fData.size()); } @@ -3910,7 +3880,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierRe, fData[0].theoryFourierRe->GetBinContent(i)); } fDataAvg.theoryFourierRe->SetBinContent(i, dval/fData.size()); } @@ -3924,7 +3894,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierIm, fData[0].theoryFourierIm->GetBinContent(i)); } fDataAvg.theoryFourierIm->SetBinContent(i, dval/fData.size()); } @@ -3938,7 +3908,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierPwr, fData[0].theoryFourierPwr->GetBinContent(i)); } fDataAvg.theoryFourierPwr->SetBinContent(i, dval/fData.size()); } @@ -3952,7 +3922,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierPhase, fData[0].theoryFourierPhase->GetBinContent(i)); } fDataAvg.theoryFourierPhase->SetBinContent(i, dval/fData.size()); } @@ -3963,15 +3933,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle()); } if (fDataAvg.diff != 0) { - if (!CalcAlignment(eTime)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: diff: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j 0) { - dval += fData[j].diff->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].diff, fData[0].diff->GetBinContent(i)); } fDataAvg.diff->SetBinContent(i, dval/fData.size()); } @@ -3985,7 +3950,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierRe, fData[0].diffFourierRe->GetBinContent(i)); } fDataAvg.diffFourierRe->SetBinContent(i, dval/fData.size()); } @@ -3999,7 +3964,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierIm, fData[0].diffFourierIm->GetBinContent(i)); } fDataAvg.diffFourierIm->SetBinContent(i, dval/fData.size()); } @@ -4013,7 +3978,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierPwr, fData[0].diffFourierPwr->GetBinContent(i)); } fDataAvg.diffFourierPwr->SetBinContent(i, dval/fData.size()); } @@ -4027,7 +3992,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; jGetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierPhase, fData[0].diffFourierPhase->GetBinContent(i)); } fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size()); } @@ -6314,130 +6279,40 @@ UInt_t PMusrCanvas::GetNeededAccuracy(PMsrParamStructure param) //-------------------------------------------------------------------------- -// CalcAlignment (private) +// GetInterpolatedValue (private) //-------------------------------------------------------------------------- /** - *

Calculates the alignment index for each data set needed to average the data. + *

search for xVal in histo. If xVal is not found exactly, interpolate and + * return the interpolated y-value. * * return: - * - true for perfect alignment - * - false for approximate alignment + * - interpolated value if xVal is within histo range, 0 otherwise. * - * \param tag to distinguish time data sets from Fourier data sets. + * \param histo pointer of the histogram + * \param xVal x-value to be looked for */ -Bool_t PMusrCanvas::CalcAlignment(const EAlignTag tag) +Double_t PMusrCanvas::GetInterpolatedValue(TH1F* histo, Double_t xVal) { - Bool_t result = true; + if (histo == 0) + return 0.0; - fAlignmentOffset.clear(); + Int_t idx = histo->FindBin(xVal); - UInt_t idx=0; - Double_t dval; - if (tag == eTime) { - // first find the data vector with the lowest initial time - dval = fData[0].data->GetXaxis()->GetBinCenter(1); - for (UInt_t i=1; iGetXaxis()->GetBinCenter(1) < dval) { - idx = i; - dval = fData[i].data->GetXaxis()->GetBinCenter(1); - } - } + // make sure idx is within range + if ((idx < 1) || (idx > histo->GetNbinsX())) + return 0.0; - // next setp: find all the alignment indices - fAlignmentOffset.resize(fData.size()); - for (UInt_t i=0; iGetXaxis()->GetBinCenter(1); - Int_t j=1; - Bool_t found = false; - do { - if (fData[idx].data->GetXaxis()->GetBinCenter(j) >= dval) { - fAlignmentOffset[i] = (UInt_t)(j-1); - if (fData[idx].data->GetXaxis()->GetBinCenter(j) != dval) { - result = false; - } - found = true; - } - } while (!found && (++jGetNbinsX())); - } - } else if (tag == eTheoTime) { - // first find the data vector with the lowest initial time - dval = fData[0].theory->GetXaxis()->GetBinCenter(1); - for (UInt_t i=1; iGetXaxis()->GetBinCenter(1) < dval) { - idx = i; - dval = fData[i].theory->GetXaxis()->GetBinCenter(1); - } - } + // make sure the lower bound idx is found. This is needed since + // FindBin rounds towards the closer idx. + if (histo->GetBinCenter(idx) > xVal) + --idx; - // next setp: find all the alignment indices - fAlignmentOffset.resize(fData.size()); - for (UInt_t i=0; iGetXaxis()->GetBinCenter(1); - Int_t j=1; - Bool_t found = false; - do { - if (fData[idx].theory->GetXaxis()->GetBinCenter(j) >= dval) { - fAlignmentOffset[i] = (UInt_t)(j-1); - if (fData[idx].theory->GetXaxis()->GetBinCenter(j) != dval) { - result = false; - } - found = true; - } - } while (!found && (++jGetNbinsX())); - } - } else if (tag == eFreq) { - // first find the data vector with the lowest initial time - dval = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(1); - for (UInt_t i=1; iGetXaxis()->GetBinCenter(1) < dval) { - idx = i; - dval = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1); - } - } + Double_t x0, x1, y0, y1; + x0 = histo->GetBinCenter(idx); + x1 = histo->GetBinCenter(idx+1); + y0 = histo->GetBinContent(idx); + y1 = histo->GetBinContent(idx+1); - // next setp: find all the alignment indices - fAlignmentOffset.resize(fData.size()); - for (UInt_t i=0; iGetXaxis()->GetBinCenter(1); - Int_t j=1; - Bool_t found = false; - do { - if (fData[idx].dataFourierRe->GetXaxis()->GetBinCenter(j) >= dval) { - fAlignmentOffset[i] = (UInt_t)(j-1); - if (fData[idx].dataFourierRe->GetXaxis()->GetBinCenter(j) != dval) { - result = false; - } - found = true; - } - } while (!found && (++jGetNbinsX())); - } - } else if (tag == eTheoFreq) { - // first find the data vector with the lowest initial time - dval = fData[0].theoryFourierRe->GetXaxis()->GetBinCenter(1); - for (UInt_t i=1; iGetXaxis()->GetBinCenter(1) < dval) { - idx = i; - dval = fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1); - } - } - - // next setp: find all the alignment indices - fAlignmentOffset.resize(fData.size()); - for (UInt_t i=0; iGetXaxis()->GetBinCenter(1); - Int_t j=1; - Bool_t found = false; - do { - if (fData[idx].theoryFourierRe->GetXaxis()->GetBinCenter(j) >= dval) { - fAlignmentOffset[i] = (UInt_t)(j-1); - if (fData[idx].theoryFourierRe->GetXaxis()->GetBinCenter(j) != dval) { - result = false; - } - found = true; - } - } while (!found && (++jGetNbinsX())); - } - } - - return result; + return (y1-y0)*(xVal-x0)/(x1-x0)+y0; } + diff --git a/src/classes/PPrepFourier.cpp b/src/classes/PPrepFourier.cpp index 30750f8f..cf077229 100644 --- a/src/classes/PPrepFourier.cpp +++ b/src/classes/PPrepFourier.cpp @@ -70,14 +70,12 @@ PPrepFourier::~PPrepFourier() // SetBkgRange //-------------------------------------------------------------------------- /** - *

+ *

set the background range. * - * \param bkgRange + * \param bkgRange array with background range */ void PPrepFourier::SetBkgRange(const Int_t *bkgRange) { - cout << endl << "debug> bkgRange: " << bkgRange[0] << ", " << bkgRange[1] << endl; - int err=0; if (bkgRange[0] >= -1) { fBkgRange[0] = bkgRange[0]; @@ -118,9 +116,9 @@ void PPrepFourier::SetBkgRange(const Int_t *bkgRange) // SetPacking //-------------------------------------------------------------------------- /** - *

+ *

set the packing for the histograms. * - * \param packing + * \param packing number to be used. */ void PPrepFourier::SetPacking(const Int_t packing) { @@ -135,9 +133,10 @@ void PPrepFourier::SetPacking(const Int_t packing) // AddData //-------------------------------------------------------------------------- /** - *

+ *

add a data-set (time domain data + meta information) to the internal + * data vector. * - * \param data + * \param data set to be added */ void PPrepFourier::AddData(musrFT_data &data) { @@ -148,13 +147,10 @@ void PPrepFourier::AddData(musrFT_data &data) // DoBkgCorrection //-------------------------------------------------------------------------- /** - *

- * + *

Correct the internal data sets according to a background interval given. */ void PPrepFourier::DoBkgCorrection() { - cout << endl << "debug> DoBkgCorrection ..."; - // make sure fData are already present, and if not create the necessary data sets if (fData.size() != fRawData.size()) { InitData(); @@ -162,7 +158,6 @@ void PPrepFourier::DoBkgCorrection() // if no bkg-range is given, nothing needs to be done if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) { - cout << endl << "debug> no background range given ..."; return; } @@ -181,7 +176,6 @@ void PPrepFourier::DoBkgCorrection() bkg += fRawData[i].rawData[j]; } bkg /= (fBkgRange[1]-fBkgRange[0]+1); - cout << endl << "debug> histo: " << i << ", bkg=" << bkg; // correct data for (unsigned int j=0; j - * + *

Rebin (pack) the internal data. */ void PPrepFourier::DoPacking() { - cout << endl << "debug> DoPacking : pack=" << fPacking << " ..."; - // make sure fData are already present, and if not create the necessary data sets if (fData.size() != fRawData.size()) { InitData(); @@ -224,10 +215,6 @@ void PPrepFourier::DoPacking() // change the original data set with the packed one fData[i].clear(); fData[i] = tmpData; - - cout << endl << "debug> histo " << i+1 << ": packed data: "; - for (unsigned int j=0; j<15; j++) - cout << fData[i][j] << ", "; } } @@ -235,27 +222,25 @@ void PPrepFourier::DoPacking() // DoFiltering //-------------------------------------------------------------------------- /** - *

- * + *

Not implemented yet. */ void PPrepFourier::DoFiltering() { - cout << endl << "debug> DoFiltering not yet implemented ..."; - // make sure fData are already present, and if not create the necessary data sets if (fData.size() != fRawData.size()) { InitData(); } - } //-------------------------------------------------------------------------- // DoLifeTimeCorrection //-------------------------------------------------------------------------- /** - *

+ *

Try to do a muon life time correction. The idea is to estimate N0 without + * any theory. This will be OK for high fields (> couple kGauss) but not so good + * for low fields. * - * \param fudge + * \param fudge rescaling factor for the estimated N0. Should be around 1 */ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge) { @@ -294,9 +279,9 @@ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge) // GetInfo //-------------------------------------------------------------------------- /** - *

+ *

Returns the meta information of a data set. * - * \param idx + * \param idx index of the object */ TString PPrepFourier::GetInfo(const UInt_t idx) { @@ -312,8 +297,8 @@ TString PPrepFourier::GetInfo(const UInt_t idx) // GetData //-------------------------------------------------------------------------- /** - *

- * + *

Creates the requested TH1F objects and returns them. The ownership is with + * the caller. */ vector PPrepFourier::GetData() { @@ -382,7 +367,7 @@ vector PPrepFourier::GetData() *

Creates the requested TH1F object and returns it. The ownership is with * the caller. * - * \param idx + * \param idx index of the requested histogram */ TH1F *PPrepFourier::GetData(const UInt_t idx) { diff --git a/src/include/PFourier.h b/src/include/PFourier.h index acb7a5b2..787450e1 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -52,6 +52,8 @@ class PFourier virtual void Transform(UInt_t apodizationTag = 0); + virtual const char* GetDataTitle() { return fData->GetTitle(); } + virtual const Int_t GetUnitTag() { return fUnitTag; } virtual Double_t GetResolution() { return fResolution; } virtual TH1F* GetRealFourier(const Double_t scale = 1.0); virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0); diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h new file mode 100644 index 00000000..0634fb96 --- /dev/null +++ b/src/include/PFourierCanvas.h @@ -0,0 +1,161 @@ +/*************************************************************************** + * Copyright (C) 2007-2015 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 _PFOURIERCANVAS_H_ +#define _PFOURIERCANVAS_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "PMusr.h" +#include "PFourier.h" + +// Canvas menu id's +#define P_MENU_ID_FOURIER 10001 +#define P_MENU_ID_AVERAGE 10002 +#define P_MENU_ID_SAVE_DATA 10003 + +#define P_MENU_ID_FOURIER_REAL 100 +#define P_MENU_ID_FOURIER_IMAG 101 +#define P_MENU_ID_FOURIER_REAL_AND_IMAG 102 +#define P_MENU_ID_FOURIER_PWR 103 +#define P_MENU_ID_FOURIER_PHASE 104 +#define P_MENU_ID_FOURIER_PHASE_PLUS 105 +#define P_MENU_ID_FOURIER_PHASE_MINUS 106 + +#define P_MENU_ID_SAVE_ASCII 200 + +//------------------------------------------------------------------------ +/** + *

Structure holding all necessary Fourier histograms. + */ +typedef struct { + TH1F *dataFourierRe; ///< real part of the Fourier transform of the data histogram + TH1F *dataFourierIm; ///< imaginary part of the Fourier transform of the data histogram + TH1F *dataFourierPwr; ///< power spectrum of the Fourier transform of the data histogram + TH1F *dataFourierPhase; ///< phase spectrum of the Fourier transform of the data histogram +} PFourierCanvasDataSet; + +//------------------------------------------------------------------------ +/** + *

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

+ */ +class PFourierCanvas : public TObject, public TQObject +{ + public: + PFourierCanvas(); + PFourierCanvas(vector &fourier, const Char_t* title, const Bool_t showAverage, + const Int_t fourierPlotOpt, Double_t fourierXrange[2], Double_t phase, + Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, const Bool_t batch); + PFourierCanvas(vector &fourier, const Char_t* title, const Bool_t showAverage, + const Int_t fourierPlotOpt, Double_t fourierXrange[2], Double_t phase, + Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, + const PIntVector markerList, const PIntVector colorList, const Bool_t batch); + virtual ~PFourierCanvas(); + + virtual void Done(Int_t status=0); // *SIGNAL* + virtual void HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected); // SLOT + virtual void HandleMenuPopup(Int_t id); // SLOT + virtual void LastCanvasClosed(); // SLOT + + virtual void UpdateFourierPad(); + virtual void UpdateInfoPad(); + + virtual Bool_t IsValid() { return fValid; } + + virtual void SetTimeout(Int_t ival); + + virtual void SaveGraphicsAndQuit(const Char_t *fileName); + virtual void SaveDataAscii(); + + private: + Int_t fTimeout; ///< timeout after which the Done signal should be emited. If timeout <= 0, no timeout is taking place + Bool_t fBatchMode; ///< musrview in ROOT batch mode + Bool_t fValid; ///< if true, everything looks OK + Bool_t fAveragedView; ///< tag showing that the averaged view or normal view should be presented. + Int_t fCurrentPlotView; ///< tag showing what the current plot view is: real, imag, power, phase, ... + Double_t fInitialXRange[2]; ///< keeps the initial x-range + Double_t fInitialYRange[2]; ///< keeps the initial y-range + + TString fTitle; + TString fXaxisTitle; + vector fFourier; ///< keeps all the Fourier data, ownership is with the caller + PFourierCanvasDataList fFourierHistos; ///< keeps all the Fourier histos + PFourierCanvasDataSet fFourierAverage; ///< keeps the average of the Fourier histos + Double_t fCurrentFourierPhase; ///< keeps the current Fourier phase (real/imag) + TLatex *fCurrentFourierPhaseText; ///< used in Re/Im Fourier to show the current phase in the pad + + TStyle *fStyle; ///< A collection of all graphics attributes + + TTimer *fTimeoutTimer; ///< timeout timer in order to terminate if no action is taking place for too long + + PIntVector fMarkerList; ///< list of markers + PIntVector fColorList; ///< list of colors + + // canvas menu related variables + TRootCanvas *fImp; ///< ROOT native GUI version of main window with menubar and drawing area + TGMenuBar *fBar; ///< menu bar + TGPopupMenu *fPopupMain; ///< popup menu MusrFT in the main menu bar + TGPopupMenu *fPopupSave; ///< popup menu of the MusrFT/Save Data sub menu + TGPopupMenu *fPopupFourier; ///< popup menu of the MusrFT/Fourier sub menu + + // canvas related variables + TCanvas *fMainCanvas; ///< main canvas + TPaveText *fTitlePad; ///< title pad used to display a title + TPad *fFourierPad; ///< fourier pad used to display the fourier + TLegend *fInfoPad; ///< info pad used to display a legend of the data plotted + + virtual void CreateXaxisTitle(); + virtual void CreateStyle(); + virtual void InitFourierDataSets(); + virtual void InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh); + virtual void CleanupAverage(); + virtual void HandleAverage(); + + virtual void PlotFourier(); + virtual void PlotFourierPhaseValue(); + virtual void PlotAverage(); + virtual void IncrementFourierPhase(); + virtual void DecrementFourierPhase(); + + virtual Double_t GetMaximum(TH1F* histo, Double_t xmin=-1.0, Double_t xmax=-1.0); + virtual Double_t GetMinimum(TH1F* histo, Double_t xmin=-1.0, Double_t xmax=-1.0); + virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal); + + ClassDef(PFourierCanvas, 1) +}; + +#endif // _PFOURIERCANVAS_H_ diff --git a/src/include/PFourierCanvasLinkDef.h b/src/include/PFourierCanvasLinkDef.h new file mode 100644 index 00000000..8ceb2d47 --- /dev/null +++ b/src/include/PFourierCanvasLinkDef.h @@ -0,0 +1,39 @@ +/*************************************************************************** + + PMusrCanvasLinkDef.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2015 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 __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class PFourierCanvas+; + +#endif + diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 7b44260e..e8296e4c 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -278,9 +278,7 @@ class PMusrCanvas : public TObject, public TQObject PRunListCollection *fRunList; ///< data handler #endif // __MAKECINT__ - enum EAlignTag {eTime, eTheoTime, eFreq, eTheoFreq}; PMusrCanvasDataSet fDataAvg; ///< set of all averaged data to be plotted (asymmetry/single histogram) - PIntVector fAlignmentOffset; ///< holds the vector with the time/freq alignment offsets PMusrCanvasDataList fData; ///< list of all histogram data to be plotted (asymmetry/single histogram) PMusrCanvasNonMusrDataList fNonMusrData; ///< list of all error graphs to be plotted (non-muSR) @@ -333,7 +331,7 @@ class PMusrCanvas : public TObject, public TQObject virtual Bool_t IsScaleN0AndBkg(); virtual UInt_t GetNeededAccuracy(PMsrParamStructure param); - virtual Bool_t CalcAlignment(const EAlignTag tag); + virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal); ClassDef(PMusrCanvas, 1) }; diff --git a/src/include/PPrepFourier.h b/src/include/PPrepFourier.h index 735a23b7..d1066809 100644 --- a/src/include/PPrepFourier.h +++ b/src/include/PPrepFourier.h @@ -40,14 +40,15 @@ using namespace std; //---------------------------------------------------------------------------- /** - *

+ *

Data structure holding raw time domain uSR data together with some + * necessary meta information. */ typedef struct { - TString info; //< keeps all the meta information - double timeResolution; //< time resolution in (usec) - int t0; //< keep the t0 bin - Double_t timeRange[2]; //< time range to be used, given in (usec). - PDoubleVector rawData; //< a single time domain data vector + TString info; ///< keeps all the meta information + double timeResolution; ///< time resolution in (usec) + int t0; ///< keep the t0 bin + Double_t timeRange[2]; ///< time range to be used, given in (usec). + PDoubleVector rawData; ///< a single time domain data vector } musrFT_data; //---------------------------------------------------------------------------- diff --git a/src/musrFT.cpp b/src/musrFT.cpp index 3551eedb..edf84510 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -38,6 +38,8 @@ #include using namespace std; +#include +#include #include #include #include @@ -50,36 +52,39 @@ using namespace std; #include "PRunDataHandler.h" #include "PPrepFourier.h" #include "PFourier.h" +#include "PFourierCanvas.h" //---------------------------------------------------------------------------- /** - *

+ *

Structure keeping the command line options. */ typedef struct { - vector msrFln; - vector dataFln; - vector dataFileFormat; - TString graphicFormat; - TString dumpFln; - int bkg[2]; - TString fourierOpt; - TString apotization; - int fourierPower; - TString fourierUnits; - double initialPhase; - double fourierRange[2]; - double timeRange[2]; - vector histo; - bool showAverage; - vector t0; - int packing; - TString title; - double lifetimecorrection; //< is == 0.0 for NO life time correction, otherwise it holds the fudge factor + vector msrFln; ///< msr-file names to be used. + vector dataFln; ///< raw-data-file names to be used. + vector dataFileFormat; ///< file format guess + TString graphicFormat; ///< format for the graphical output dump + TString dumpFln; ///< dump file name for Fourier data output + TString msrFlnOut; ///< dump file name for msr-file generation + int bkg[2]; ///< background range + TString fourierOpt; ///< Fourier options, i.e. real, imag, power, phase + TString apodization; ///< apodization setting: none, weak, medium, strong + int fourierPower; ///< Fourier power for zero padding, i.e. 2^fourierPower points + TString fourierUnits; ///< wished Fourier units: Gauss, Tesla, MHz, Mc/s + double initialPhase; ///< inital Fourier phase for Real/Imag + double fourierRange[2]; ///< Fourier range to be plotted. Given in the choosen units. + double timeRange[2]; ///< time range used for the Fourier + vector histo; ///< selection of the histos used from at data file for Fourier + bool showAverage; ///< flag indicating if initially the Fourier average over the given histos shall be plotted. + vector t0; ///< t0 vector for the histos. If not given t0's will be estimated. + int packing; ///< packing for rebinning the time histograms before Fourier transform. + TString title; ///< title to be shown for the Fourier plot. + double lifetimecorrection; ///< is == 0.0 for NO life time correction, otherwise it holds the fudge factor + Int_t timeout; ///< timeout in (sec) after which musrFT will terminate. if <= 0, no automatic termination will take place. } musrFT_startup_param; //------------------------------------------------------------------------- /** - *

+ *

prints the musrFT usage. */ void musrFT_syntax() { @@ -104,9 +109,9 @@ void musrFT_syntax() cout << endl << " --filter : filter and filter-specific-information -- ***TO BE WRITTEN YET***."; cout << endl << " -b, --background : background interval used to estimate the backround to be"; cout << endl << " subtracted before the Fourier transform. , to be given in bins."; - cout << endl << " -fo, --fourier-option : can be 'real', 'power', 'imag', 'real+imag', of 'phase'."; + cout << endl << " -fo, --fourier-option : can be 'real', 'imag', 'real+imag', 'power', or 'phase'."; cout << endl << " If this is not defined (neither on the command line nor in the musrFT_startup.xml),"; - cout << endl << " default will be 'real'."; + cout << endl << " default will be 'power'."; cout << endl << " -apod, --apodization : can be either 'none', 'weak', 'medium', 'strong'."; cout << endl << " Default will be 'none'."; cout << endl << " -fp, --fourier-power : being the Fourier power, i.e. 2^ used for zero padding."; @@ -115,7 +120,7 @@ void musrFT_syntax() cout << endl << " One may choose between the fields (Gauss) or (Tesla), the frequency (MHz),"; cout << endl << " and the angular-frequency domain (Mc/s)."; cout << endl << " Default will be 'MHz'."; - cout << endl << " -ph, --phase : defines the inital phase . This only is of concern for 'real',"; + cout << endl << " -ph, --phase : defines the initial phase . This only is of concern for 'real',"; cout << endl << " '', and 'real+imag'."; cout << endl << " Default will be 0.0."; cout << endl << " -fr, --fourier-range : Fourier range. , are interpreted in the units given."; @@ -134,12 +139,14 @@ void musrFT_syntax() cout << endl << " to all histos."; cout << endl << " Example: musrFT -df lem15_his_01234.root -fo real --t0 2750 --histo 1 3"; cout << endl << " -pa, --packing : if (an integer), the time domain data will first be packed/rebinned by ."; - cout << endl << " -t, --title : give a global title for the plot."; + cout << endl << " --title <title> : give a global title for the plot."; cout << endl << " --create-msr-file <fln> : creates a msr-file based on the command line options"; cout << endl << " provided. This will help on the way to a full fitting model."; cout << endl << " ***TO BE WRITTEN YET.***"; cout << endl << " -lc, --lifetimecorrection <fudge>: try to eliminate muon life time decay. Only makes sense for low"; cout << endl << " transverse fields. <fudge> is a tweaking factor and should be kept around 1.0."; + cout << endl << " --timeout <timeout> : <timeout> given in seconds after which musrFT terminates."; + cout << endl << " If <timeout> <= 0, no timeout will take place. Default <timeout> is 3600."; cout << endl << endl; } @@ -147,16 +154,17 @@ void musrFT_syntax() /** * <p>initialize startup parameters. * - * \param startupParam + * \param startupParam command line options */ void musrFT_init(musrFT_startup_param &startupParam) { startupParam.graphicFormat = TString(""); startupParam.dumpFln = TString(""); + startupParam.msrFlnOut = TString(""); startupParam.bkg[0] = -1; startupParam.bkg[1] = -1; - startupParam.fourierOpt = TString("real"); - startupParam.apotization = TString("none"); + startupParam.fourierOpt = TString("??"); + startupParam.apodization = TString("none"); startupParam.fourierPower = -1; startupParam.fourierUnits = TString("??"); startupParam.initialPhase = 0.0; @@ -168,6 +176,7 @@ void musrFT_init(musrFT_startup_param &startupParam) startupParam.packing = 1; startupParam.title = TString(""); startupParam.lifetimecorrection = 0.0; + startupParam.timeout = 3600; } //------------------------------------------------------------------------- @@ -273,9 +282,9 @@ bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &s * * <b>return:</b> 0 if everything is OK, 1 for --version or --help, 2 for an error. * - * \param argc - * \param argv - * \param startupParam + * \param argc number of command line arguments + * \param argv command line argument array + * \param startupParam command line data structure */ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupParam) { @@ -356,7 +365,7 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa cerr << endl << ">> musrFT **ERROR** found option --apodization with unrecognized argument '" << topt << "'." << endl; return 2; } - startupParam.apotization = topt; + startupParam.apodization = topt; i++; } else if (tstr.BeginsWith("-fp") || tstr.BeginsWith("--fourier-power")) { if (i+1 >= argc) { // something is wrong since there needs to be two arguments here @@ -449,7 +458,7 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa cerr << endl << ">> musrFT **ERROR** found option --t0 without argument!" << endl; return 2; } - } else if (tstr.BeginsWith("-t ") || tstr.BeginsWith("--title")) { + } else if (tstr.BeginsWith("--title")) { if (i+1 >= argc) { // something is wrong since there needs to be an argument here cerr << endl << ">> musrFT **ERROR** found option --title without argument!" << endl; return 2; @@ -468,6 +477,13 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa return 2; } startupParam.packing = pack.Atoi(); + } else if (tstr.BeginsWith("--create-msr-file")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --create-msr-file without argument!" << endl; + return 2; + } + ++i; + startupParam.msrFlnOut = TString(argv[i]); } else if (tstr.BeginsWith("-lc") || tstr.BeginsWith("--lifetimecorrection")) { if (i+1 >= argc) { // something is wrong since there needs to be an argument here cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection without argument!" << endl; @@ -476,10 +492,22 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa ++i; TString fudge(argv[i]); if (!fudge.IsFloat()) { - cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection with a fudge which is not an double '" << fudge << "'." << endl; + cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection with a fudge which is not a double '" << fudge << "'." << endl; return 2; } startupParam.lifetimecorrection = fudge.Atof(); + } else if (tstr.BeginsWith("--timeout")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --timeout without argument!" << endl; + return 2; + } + ++i; + TString tt(argv[i]); + if (!tt.IsDigit()) { + cerr << endl << ">> musrFT **ERROR** found option --timeout with a <timeout> which is not an integer '" << tt << "'." << endl; + return 2; + } + startupParam.timeout = tt.Atoi(); } else if (tstr.BeginsWith("-df") || tstr.BeginsWith("--data-file")) { while (++i < argc) { if (argv[i][0] == '-') { @@ -541,7 +569,11 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa //---------------------------------------------------------------------------------------- /** + * <p>Collects the meta information form the raw-data-file. * + * \param fln file name of the raw-data-file + * \param rawRunData raw-data-file object + * \param metaInfo return string which will contain the meta information. */ void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &metaInfo) { @@ -581,6 +613,13 @@ void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &met } //------------------------------------------------------------------------- +/** + * <p>Estimates the t0's of the raw-data-files. It simply is looking for the + * maximum of the raw-data (assuming a prompt peak). This will fail for LEM + * and ISIS data for sure. + * + * \param rd raw-data-file collection (see PPrepFourier.h) + */ void musrFT_estimateT0(musrFT_data &rd) { cout << endl << ">> musrFT **WARNING** try to estimate t0 from maximum in the data set"; @@ -600,12 +639,12 @@ void musrFT_estimateT0(musrFT_data &rd) } //------------------------------------------------------------------------- -void musrFT_cleanup(TH1F *h) /** - * <p> + * <p> deletes a histogram. * - * \param h histogram to be deleted + * \param h point to a ROOT histogram object */ +void musrFT_cleanup(TH1F *h) { if (h) { delete h; @@ -615,9 +654,12 @@ void musrFT_cleanup(TH1F *h) //------------------------------------------------------------------------- /** - * <p> + * <p>Dump the Fourier transformed data into an ascii file. * - * \param data + * \param fln dump file name + * \param fourierData collection of all the Fourier transformed data. + * \param start starting point from where the data shall be written to file. + * \param end ending point up to where the data shall be written to file. */ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, double end) { @@ -691,14 +733,15 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d return 0; } -//--------------------------------------------------- +//------------------------------------------------------------------------- /** - * <p> + * <p>Groups the histograms before Fourier transform. This is used to group + * detectors. * - * \param runDataHandler - * \param global - * \param run - * \param rd + * \param runDataHandler raw-run-data object containing the data + * \param global pointer to the GLOBAL block of the msr-file + * \param run reference to the relevant RUN block of the msr-file + * \param rd data collection which will hold the grouped histograms. */ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, PMsrRunBlock &run, musrFT_data &rd) { @@ -779,7 +822,115 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, return 0; } -//--------------------------------------------------- +//------------------------------------------------------------------------- +/** + * <p>Dumps an msr-file according to the given command line settings. This + * is meant to generate an initial msr-file for a given data-file. This + * routine is 'stupid' in the sense that it knows nothing about the data-files. + * Hence when feeding it with senseless command line settings, the resulting + * msr-file fed back to musrFT might do funny things! + * + * \param param command line options + */ +void musrFT_dumpMsrFile(musrFT_startup_param ¶m) +{ + ofstream fout(param.msrFlnOut.Data(), ofstream::out); + + // write title + if (param.title.Length() == 0) { // create title if not given + if (param.dataFln.size() != 0) { + param.title = param.dataFln[0]; + } else { + param.title = param.msrFlnOut; + } + } + fout << param.title << endl; + fout << "###############################################################" << endl; + + // write GLOBAL block + fout << "GLOBAL" << endl; + fout << "fittype 0 (single histogram fit)" << endl; + if (param.t0.size() == 1) { // only a single t0 value given, hence assume it is valid for ALL histos + fout << "t0 " << param.t0[0] << endl; + } + if ((param.timeRange[0] != -1.0) && (param.timeRange[1] != -1.0)) { + fout << "fit " << param.timeRange[0] << " " << param.timeRange[1] << endl; + } + fout << "packing " << param.packing << endl; + fout << endl; + fout << "###############################################################" << endl; + + // write RUN block + // get extension of the data file + TString fileFormat("MUSR-ROOT"); + for (unsigned int i=0; i<param.dataFln.size(); i++) { + if (param.dataFileFormat[i].BeginsWith("PsiBin")) + fileFormat = TString("PSI-MDU"); + else if (param.dataFileFormat[i].BeginsWith("NeXus")) + fileFormat = TString("NEXUS"); + else if (param.dataFileFormat[i].BeginsWith("Mud")) + fileFormat = TString("MUD"); + for (unsigned int j=0; j<param.histo.size(); j++) { + fout << "RUN " << param.dataFln[i] << " BXXX IXX " << fileFormat << " (name beamline institute data-file-format)" << endl; + fout << "forward " << param.histo[j] << endl; + if ((param.t0.size() > 1) && (j < param.t0.size())) { + fout << "t0 " << param.t0[j] << endl; + } + if ((param.bkg[0] > -1) && (param.bkg[1] > -1)) + fout << "background " << param.bkg[0] << " " << param.bkg[1] << endl; + fout << "#--------------------------------------------------------------" << endl; + } + } + fout << endl; + fout << "###############################################################" << endl; + + // write PLOT block + fout << "PLOT 0 (single histo plot)" << endl; + if (param.histo.size() == 0) { + fout << "runs 1" << endl; + } else { + fout << "runs "; + for (unsigned int i=0; i<param.histo.size(); i++) + fout << i+1 << " "; + fout << endl; + } + if ((param.timeRange[0] == -1.0) && (param.timeRange[1] == -1.0)) { + fout << "range 0 10" << endl; + } else { + fout << "range " << param.timeRange[0] << " " << param.timeRange[1] << endl; + } + fout << endl; + fout << "###############################################################" << endl; + + // write FOURIER block + fout << "FOURIER" << endl; + if (param.fourierUnits.BeginsWith("??")) { // Fourier units not given, hence choose MHz + fout << "units MHz # units either 'Gauss', 'MHz', or 'Mc/s'" << endl; + } else { + fout << "units " << param.fourierUnits << " # units either 'Gauss', 'MHz', or 'Mc/s'" << endl; + } + if (param.fourierOpt.BeginsWith("??")) { // Fourier plot option not given, hence choose POWER + fout << "plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE" << endl; + } else { + fout << "plot " << param.fourierOpt << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE" << endl; + } + if (param.fourierPower > 1) { + fout << "fourier_power " << param.fourierPower << endl; + } + fout << "apodization " << param.apodization << " # NONE, WEAK, MEDIUM, STRONG" << endl; + if ((param.fourierRange[0] > -1.0) && (param.fourierRange[1] > -1.0)) { + fout << "range " << param.fourierRange[0] << " " << param.fourierRange[1] << endl; + } + + fout.close(); +} + +//------------------------------------------------------------------------- +/** + * <p>Gets time a time stamp in msec. Used to measure the calculation time. + * + * <b>return:</b> time stamp with msec resolution. + */ double millitime() { struct timeval now; @@ -790,15 +941,20 @@ double millitime() //------------------------------------------------------------------------- /** - * <p> + * <p>musrFT is used to do a Fourier transform of uSR data without any fitting. + * It directly Fourier transforms the raw histogram data (exception see --lifetimecorrection), + * and hence only will give staisfactory results for applied fields of larger a + * couple of kGauss. It is meant to be used to get a feeling what time-domain + * model will be appropriate. It is NOT meant for ANY quantitative analysis! * - * \param argc - * \param argv + * \param argc number of command line arguments + * \param argv command line argument array */ int main(int argc, char *argv[]) { Int_t unitTag = FOURIER_UNIT_NOT_GIVEN; Int_t apodTag = F_APODIZATION_NONE; + Int_t fourierPlotTag = FOURIER_PLOT_NOT_GIVEN; // only program name alone if (argc == 1) { @@ -821,6 +977,12 @@ int main(int argc, char *argv[]) return retVal; } + // dump msr-file + if (startupParam.msrFlnOut.Length() > 0) { + musrFT_dumpMsrFile(startupParam); + return PMUSR_SUCCESS; + } + // read startup file char startup_path_name[128]; PStartupOptions startup_options; @@ -903,7 +1065,6 @@ int main(int argc, char *argv[]) // load data-file(s) provided directly for (unsigned int i=msrHandler.size(); i<msrHandler.size()+startupParam.dataFln.size(); i++) { -// cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()]; // create run data handler if (startupHandler) runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList()); @@ -956,6 +1117,9 @@ int main(int argc, char *argv[]) unsigned int idx=0; // get meta info, time resolution, time range, raw data sets if (i < msrHandler.size()) { // obtain info from msr-files + // keep title if not overwritten by the command line + if (startupParam.title.Length() == 0) + startupParam.title = *(msrHandler[0]->GetMsrTitle()); // keep PLOT block info PMsrPlotList *plot = msrHandler[i]->GetMsrPlotList(); if (plot == 0) { @@ -993,7 +1157,11 @@ int main(int argc, char *argv[]) // get range if ((startupParam.fourierRange[0] == -1) && (startupParam.fourierRange[1] == -1)) { // no Fourier range given from the command line startupParam.fourierRange[0] = fourierBlock->fPlotRange[0]; - startupParam.fourierRange[1] = fourierBlock->fPlotRange[1]; + startupParam.fourierRange[1] = fourierBlock->fPlotRange[1]; + } + // get Fourier plot option, i.e. real, imag, power, phase + if (startupParam.fourierOpt.BeginsWith("??")) { // only do something if not overwritten by the command line + fourierPlotTag = fourierBlock->fPlotTag; } } } @@ -1030,8 +1198,10 @@ int main(int argc, char *argv[]) rd.timeRange[0] = startupParam.timeRange[0]; rd.timeRange[1] = startupParam.timeRange[1]; } else { - rd.timeRange[0] = plot->at(0).fTmin[0]; - rd.timeRange[1] = plot->at(0).fTmax[0]; + if (plot->at(0).fTmin.size() > 0) { + rd.timeRange[0] = plot->at(0).fTmin[0]; + rd.timeRange[1] = plot->at(0).fTmax[0]; + } } // handle data set(s) @@ -1089,6 +1259,22 @@ int main(int argc, char *argv[]) } } + // make sure Fourier plot tag is set + if (fourierPlotTag == FOURIER_PLOT_NOT_GIVEN) { + if (!startupParam.fourierOpt.CompareTo("real", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_REAL; + else if (!startupParam.fourierOpt.CompareTo("imag", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_IMAG; + else if (!startupParam.fourierOpt.CompareTo("real+imag", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_REAL_AND_IMAG; + else if (!startupParam.fourierOpt.CompareTo("power", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_POWER; + else if (!startupParam.fourierOpt.CompareTo("phase", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_PHASE; + else + fourierPlotTag = FOURIER_PLOT_POWER; + } + // calculate background levels and subtract them from the data data.DoBkgCorrection(); @@ -1124,15 +1310,13 @@ int main(int argc, char *argv[]) } // Fourier transform data - if (startupParam.apotization.BeginsWith("weak", TString::kIgnoreCase)) + if (startupParam.apodization.BeginsWith("weak", TString::kIgnoreCase)) apodTag = F_APODIZATION_WEAK; - else if (startupParam.apotization.BeginsWith("medium", TString::kIgnoreCase)) + else if (startupParam.apodization.BeginsWith("medium", TString::kIgnoreCase)) apodTag = F_APODIZATION_MEDIUM; - else if (startupParam.apotization.BeginsWith("strong", TString::kIgnoreCase)) + else if (startupParam.apodization.BeginsWith("strong", TString::kIgnoreCase)) apodTag = F_APODIZATION_STRONG; - cout << endl << "debug> apodTag = " << apodTag << endl; - double start = millitime(); for (unsigned int i=0; i<fourier.size(); i++) { fourier[i]->Transform(apodTag); @@ -1140,19 +1324,83 @@ int main(int argc, char *argv[]) double end = millitime(); cout << endl << "debug> after FFT. calculation time: " << (end-start)/1.0e3 << " (sec)." << endl; + PFourierCanvas *fourierCanvas = 0; + // if Fourier dumped if whished do it now if (startupParam.dumpFln.Length() > 0) { musrFT_dumpData(startupParam.dumpFln, fourier, startupParam.fourierRange[0], startupParam.fourierRange[1]); - } else { + } else { // do Canvas - // if Fourier graphical export is whished, switch to batch mode + // if Fourier graphical export is whished, switch to batch mode + Bool_t batch = false; + if (startupParam.graphicFormat.Length() != 0) { + batch = true; + argv[argc] = (char*)malloc(16*sizeof(char)); + strcpy(argv[argc], "-b"); + argc++; + } - // plot the Fourier transform + // plot the Fourier transform + TApplication app("App", &argc, argv); + if (startupHandler) { + fourierCanvas = new PFourierCanvas(fourier, startupParam.title.Data(), + startupParam.showAverage, fourierPlotTag, + startupParam.fourierRange, startupParam.initialPhase, + 10, 10, 800, 800, + startupHandler->GetMarkerList(), + startupHandler->GetColorList(), + batch); + } else { + fourierCanvas = new PFourierCanvas(fourier, startupParam.title.Data(), + startupParam.showAverage, fourierPlotTag, + startupParam.fourierRange, startupParam.initialPhase, + 10, 10, 800, 800, + batch); + } + + fourierCanvas->UpdateFourierPad(); + fourierCanvas->UpdateInfoPad(); + + Bool_t ok = true; + if (!fourierCanvas->IsValid()) { + cerr << endl << ">> musrFT **SEVERE ERROR** Couldn't invoke all necessary objects, will quit."; + cerr << endl; + ok = false; + } else { + // connect signal/slot + TQObject::Connect("TCanvas", "Closed()", "PFourierCanvas", fourierCanvas, "LastCanvasClosed()"); + + fourierCanvas->SetTimeout(startupParam.timeout); + + fourierCanvas->Connect("Done(Int_t)", "TApplication", &app, "Terminate(Int_t)"); + + if (startupParam.graphicFormat.Length() != 0) { + TString fileName(""); + // create output filename based on the msr- or raw-data-filename + if (startupParam.dataFln.size() > 0) { + fileName = startupParam.dataFln[0]; + } + if (startupParam.msrFln.size() > 0) { + fileName = startupParam.msrFln[0]; + } + Ssiz_t idx = fileName.Last('.'); + fileName.Remove(idx, fileName.Length()); + fileName += "."; + fileName += startupParam.graphicFormat; + fourierCanvas->SaveGraphicsAndQuit(fileName.Data()); + } + } + // check that everything is ok + if (ok) + app.Run(true); // true needed that Run will return after quit so that cleanup works } // cleanup + if (fourierCanvas) + delete fourierCanvas; + if (startupHandler) delete startupHandler; diff --git a/src/musrview.cpp b/src/musrview.cpp index 6e035945..17d3d978 100644 --- a/src/musrview.cpp +++ b/src/musrview.cpp @@ -340,9 +340,9 @@ int main(int argc, char *argv[]) } if (asciiOutput) { - // save data in batch mode + // save data in batch mode musrCanvas->SaveDataAscii(); - musrCanvas->Done(0); + musrCanvas->Done(0); } // keep musrCanvas objects @@ -360,7 +360,6 @@ int main(int argc, char *argv[]) sprintf(canvasName, "fMainCanvas%d", i); if (gROOT->GetListOfCanvases()->FindObject(canvasName) != 0) { canvasVector[i]->~PMusrCanvas(); - } else { } } canvasVector.empty();