diff --git a/ChangeLog b/ChangeLog index ee40e460..2bd3a9c1 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,6 +4,10 @@ changes since 0.13.0 =================================== +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 2015-01-17 adding a branch for ROOT 6.x. This needs some minor adaptations due to the new rootcint/rootclang and the stricter c++11. NEW 2014-12-18 first implementation of a GLOBAL block which allows to shorten diff --git a/README b/README index ceff445b..fd52af9c 100644 --- a/README +++ b/README @@ -1,33 +1,32 @@ -# musrfit - muSR data analysis package # - -### Contents ### - -This is a data analysis package to analyze time differential muSR data. -Currently it allows the following things: - -* setting up most commonly used fitting functions for muSR -* fitting data, including global fits -* showing the fit results and the residuals -* showing the Fourier transform of the data -* extracting easily the fitting parameters to be used in other programs (gnuplot, qtiplot/origin, ...) -* allows to generate fitting input files for follow-up runs -* allows to generate global fitting input files based on a single run template -* allows to implement more sophisticated user functions - (e.g. GL vortex lattice, Meissner screening including low-energy muon stopping profiles) - -### Currently supported platforms: ### - -* Linux -* Mac OS X -* Windows - not really, only for the very brave ones - -### Documentation #### - -For a more exhaustive user documentation see: - - http://lmu.web.psi.ch/musrfit/user/MUSR/WebHome.html - -### Contact ### - - - +# musrfit - muSR data analysis package # + +### Contents ### + +This is a data analysis package to analyze time differential muSR and beta-NMR data. +Currently it allows the following things: + +* setting up most commonly used fitting functions for muSR and beta-NMR +* fitting data, including global fits +* showing the fit results and the residuals +* showing the Fourier transform of the data +* extracting easily the fitting parameters to be used in other programs (gnuplot, qtiplot/origin, ...) +* allows to generate fitting input files for follow-up runs +* allows to generate global fitting input files based on a single run template +* allows to implement more sophisticated user functions + (e.g. GL vortex lattice, Meissner screening including low-energy muon stopping profiles) + +### Currently supported platforms: ### + +* Linux +* Mac OS X +* Windows - not really, only for the very brave ones + +### Documentation #### + +For a more exhaustive user documentation see: + + http://lmu.web.psi.ch/musrfit/user/MUSR/WebHome.html + +### Contact ### + + \ No newline at end of file diff --git a/src/Makefile.am b/src/Makefile.am index 5eb0e3da..3d0d9c75 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -34,13 +34,14 @@ SUBDIRS += $(EDITORDIR) EXTRA_DIST = $(EDITORDIR)/Makefile endif -bin_PROGRAMS = musrfit musrview musrt0 msr2msr msr2data any2many +bin_PROGRAMS = musrfit musrview musrt0 musrFT msr2msr msr2data any2many bin_PROGRAMS += write_musrRoot_runHeader musrRootValidation bin_PROGRAMS += dump_header musrfit_SOURCES = musrfit.cpp musrview_SOURCES = musrview.cpp musrt0_SOURCES = musrt0.cpp +musrFT_SOURCES = musrFT.cpp msr2msr_SOURCES = msr2msr.cpp msr2data_SOURCES = msr2data.cpp any2many_SOURCES = any2many.cpp diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am index c026714d..a2d7e0cf 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 \ @@ -12,6 +13,7 @@ h_sources = \ ../include/PMusrCanvas.h \ ../include/PMusr.h \ ../include/PMusrT0.h \ + ../include/PPrepFourier.h \ ../include/PRunAsymmetry.h \ ../include/PRunBase.h \ ../include/PRunDataHandler.h \ @@ -26,6 +28,7 @@ h_sources_userFcn = \ ../include/PUserFcnBase.h h_linkdef = \ + ../include/PFourierCanvasLinkDef.h \ ../include/PMusrCanvasLinkDef.h \ ../include/PMusrT0LinkDef.h \ ../include/PStartupHandlerLinkDef.h @@ -34,6 +37,7 @@ h_linkdef_userFcn = \ ../include/PUserFcnBaseLinkDef.h dict_h_sources = \ + PFourierCanvasDict.h \ PMusrCanvasDict.h \ PMusrT0Dict.h \ PStartupHandlerDict.h @@ -45,6 +49,7 @@ cpp_sources = \ PFitter.cpp \ PFitterFcn.cpp \ PFourier.cpp \ + PFourierCanvas.cpp \ PFunction.cpp \ PFunctionHandler.cpp \ PMsr2Data.cpp \ @@ -52,6 +57,7 @@ cpp_sources = \ PMusrCanvas.cpp \ PMusr.cpp \ PMusrT0.cpp \ + PPrepFourier.cpp \ PRunAsymmetry.cpp \ PRunBase.cpp \ PRunDataHandler.cpp \ @@ -66,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/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 0fbba55a..9a91f134 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -269,7 +269,7 @@ Int_t PMsrHandler::ReadMsrFile() } // fill parameter-in-use vector - if (result == PMUSR_SUCCESS) + if ((result == PMUSR_SUCCESS) && !fFourierOnly) FillParameterInUse(theory, functions, run); // check that each run fulfills the minimum requirements @@ -278,7 +278,7 @@ Int_t PMsrHandler::ReadMsrFile() result = PMUSR_MSR_SYNTAX_ERROR; // check that parameter names are unique - if (result == PMUSR_SUCCESS) { + if ((result == PMUSR_SUCCESS) && !fFourierOnly) { UInt_t parX, parY; if (!CheckUniquenessOfParamNames(parX, parY)) { cerr << endl << ">> PMsrHandler::ReadMsrFile: **SEVERE ERROR** parameter name " << fParam[parX].fName.Data() << " is identical for parameter no " << fParam[parX].fNo << " and " << fParam[parY].fNo << "!"; @@ -290,7 +290,7 @@ Int_t PMsrHandler::ReadMsrFile() // check that if maps are present in the theory- and/or function-block, // that there are really present in the run block - if (result == PMUSR_SUCCESS) + if ((result == PMUSR_SUCCESS) && !fFourierOnly) if (!CheckMaps()) result = PMUSR_MSR_SYNTAX_ERROR; @@ -2387,6 +2387,10 @@ Int_t PMsrHandler::ParameterInUse(UInt_t paramNo) */ Bool_t PMsrHandler::HandleFitParameterEntry(PMsrLines &lines) { + // If msr-file is used for musrFT only, nothing needs to be done here. + if (fFourierOnly) + return true; + PMsrParamStructure param; Bool_t error = false; @@ -2599,6 +2603,10 @@ Bool_t PMsrHandler::HandleFitParameterEntry(PMsrLines &lines) */ Bool_t PMsrHandler::HandleTheoryEntry(PMsrLines &lines) { + // If msr-file is used for musrFT only, nothing needs to be done here. + if (fFourierOnly) + return true; + // store the theory lines fTheory = lines; @@ -2619,6 +2627,10 @@ Bool_t PMsrHandler::HandleTheoryEntry(PMsrLines &lines) */ Bool_t PMsrHandler::HandleFunctionsEntry(PMsrLines &lines) { + // If msr-file is used for musrFT only, nothing needs to be done here. + if (fFourierOnly) + return true; + // store the functions lines fFunctions = lines; @@ -3157,11 +3169,13 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines) } } // check map entries, i.e. if the map values are within parameter bounds - for (UInt_t i=0; isize(); i++) { - if ((param.GetMap(i) < 0) || (param.GetMap(i) > (Int_t) fParam.size())) { - cerr << endl << ">> PMsrHandler::HandleRunEntry: **SEVERE ERROR** map value " << param.GetMap(i) << " in line " << iter->fLineNo << " is out of range!"; - error = true; - break; + if (!fFourierOnly) { + for (UInt_t i=0; isize(); i++) { + if ((param.GetMap(i) < 0) || (param.GetMap(i) > (Int_t) fParam.size())) { + cerr << endl << ">> PMsrHandler::HandleRunEntry: **SEVERE ERROR** map value " << param.GetMap(i) << " in line " << iter->fLineNo << " is out of range!"; + error = true; + break; + } } } } @@ -3515,9 +3529,13 @@ Bool_t PMsrHandler::FilterNumber(TString str, const Char_t *filter, Int_t offset */ Bool_t PMsrHandler::HandleCommandsEntry(PMsrLines &lines) { + // If msr-file is used for musrFT only, nothing needs to be done here. + if (fFourierOnly) + return true; + PMsrLines::iterator iter; - if (lines.empty() && !fFourierOnly) { + if (lines.empty()) { cerr << endl << ">> PMsrHandler::HandleCommandsEntry(): **WARNING**: There is no COMMAND block! Do you really want this?"; cerr << endl; } @@ -4352,7 +4370,11 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) */ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines) { - if (lines.empty() && !fFourierOnly) { + // If msr-file is used for musrFT only, nothing needs to be done here. + if (fFourierOnly) + return true; + + if (lines.empty()) { cerr << endl << ">> PMsrHandler::HandleStatisticEntry: **WARNING** There is no STATISTIC block! Do you really want this?"; cerr << endl; fStatistic.fValid = false; @@ -5042,7 +5064,7 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity() } } // check fit range - if (!fRuns[i].IsFitRangeInBin()) { // fit range given as times in usec (RUN block) + if (!fRuns[i].IsFitRangeInBin() && !fFourierOnly) { // fit range given as times in usec (RUN block) if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in RUN block if (!fGlobal.IsFitRangeInBin()) { // fit range given as times in usec (GLOBAL block) if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in GLOBAL block diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 17e34378..4fcf2bd5 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1361,6 +1361,922 @@ void PMusrCanvas::SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat) Done(0); } +//-------------------------------------------------------------------------- +// SaveDataAscii +//-------------------------------------------------------------------------- +/** + *

Saves the currently seen data (data, difference, Fourier spectra, ...) in ascii column format. + */ +void PMusrCanvas::SaveDataAscii() +{ + // collect relevant data + PMusrCanvasAsciiDump dump; + PMusrCanvasAsciiDumpVector dumpVector; + + Int_t xminBin; + Int_t xmaxBin; + Double_t xmin; + Double_t xmax; + Double_t time, freq; + Double_t xval, yval; + + switch (fPlotType) { + case MSR_PLOT_SINGLE_HISTO: + case MSR_PLOT_ASYM: + case MSR_PLOT_MU_MINUS: + if (fDifferenceView) { // difference view plot + switch (fCurrentPlotView) { + case PV_DATA: + // get current x-range + xminBin = fHistoFrame->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fHistoFrame->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fHistoFrame->GetXaxis()->GetBinCenter(xminBin); + xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get time + time = fData[i].diff->GetBinCenter(j); + // check if time is in the current range + if ((time >= xmin) && (time <= xmax)) { + dump.dataX.push_back(time); + dump.data.push_back(fData[i].diff->GetBinContent(j)); + dump.dataErr.push_back(fData[i].diff->GetBinError(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_REAL: + // get current x-range + xminBin = fData[0].diffFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].diffFourierRe->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_IMAG: + // get current x-range + xminBin = fData[0].diffFourierIm->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierIm->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].diffFourierIm->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_REAL_AND_IMAG: + // get current x-range + xminBin = fData[0].diffFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].diffFourierRe->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j)); + } + } + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].diffFourierIm->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_PWR: + // get current x-range + xminBin = fData[0].diffFourierPwr->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierPwr->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].diffFourierPwr->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].diffFourierPwr->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_PHASE: + // get current x-range + xminBin = fData[0].diffFourierPhase->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierPhase->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].diffFourierPhase->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].diffFourierPhase->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + default: + break; + } + } else { // not a difference view plot + switch (fCurrentPlotView) { + case PV_DATA: + // get current x-range + xminBin = fHistoFrame->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fHistoFrame->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fHistoFrame->GetXaxis()->GetBinCenter(xminBin); + xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get time + time = fData[i].data->GetBinCenter(j); + // check if time is in the current range + if ((time >= xmin) && (time <= xmax)) { + dump.dataX.push_back(time); + dump.data.push_back(fData[i].data->GetBinContent(j)); + dump.dataErr.push_back(fData[i].data->GetBinError(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get time + time = fData[i].theory->GetBinCenter(j); + // check if time is in the current range + if ((time >= xmin) && (time <= xmax)) { + dump.theoryX.push_back(time); + dump.theory.push_back(fData[i].theory->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + + break; + case PV_FOURIER_REAL: + // get current x-range + xminBin = fData[0].dataFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].dataFourierRe->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].theoryFourierRe->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.theoryX.push_back(freq); + dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_IMAG: + // get current x-range + xminBin = fData[0].dataFourierIm->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierIm->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].dataFourierIm->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].theoryFourierIm->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.theoryX.push_back(freq); + dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_REAL_AND_IMAG: + // get current x-range + xminBin = fData[0].dataFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].dataFourierRe->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].theoryFourierRe->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.theoryX.push_back(freq); + dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + + //----------------------------- + // Im + //----------------------------- + // clean up dump + dump.dataX.clear(); + dump.data.clear(); + dump.dataErr.clear(); + dump.theoryX.clear(); + dump.theory.clear(); + + // go through all data bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].dataFourierIm->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].theoryFourierIm->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.theoryX.push_back(freq); + dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + + } + break; + case PV_FOURIER_PWR: + // get current x-range + xminBin = fData[0].dataFourierPwr->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierPwr->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].dataFourierPwr->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].dataFourierPwr->GetBinContent(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].theoryFourierPwr->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.theoryX.push_back(freq); + dump.theory.push_back(fData[i].theoryFourierPwr->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + case PV_FOURIER_PHASE: + // get current x-range + xminBin = fData[0].dataFourierPhase->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierPhase->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetNbinsX(); j++) { + // get frequency + freq = fData[i].dataFourierPhase->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.dataX.push_back(freq); + dump.data.push_back(fData[i].dataFourierPhase->GetBinContent(j)); + } + } + + // go through all theory bins + for (Int_t j=1; jGetNbinsX(); j++) { + // get frequency + freq = fData[i].theoryFourierPhase->GetBinCenter(j); + // check if time is in the current range + if ((freq >= xmin) && (freq <= xmax)) { + dump.theoryX.push_back(freq); + dump.theory.push_back(fData[i].theoryFourierPhase->GetBinContent(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + break; + default: + break; + } + } + break; + case MSR_PLOT_NON_MUSR: + if (fDifferenceView) { // difference view plot + switch (fCurrentPlotView) { + case PV_DATA: + // get current x-range + xminBin = fMultiGraphData->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fMultiGraphData->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fMultiGraphData->GetXaxis()->GetBinCenter(xminBin); + xmax = fMultiGraphData->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetN(); j++) { + // get x and y value + fNonMusrData[i].diff->GetPoint(j,xval,yval); + // check if time is in the current range + if ((xval >= xmin) && (xval <= xmax)) { + dump.dataX.push_back(xval); + dump.data.push_back(yval); + dump.dataErr.push_back(fNonMusrData[i].diff->GetErrorY(j)); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + + break; + case PV_FOURIER_REAL: + break; + case PV_FOURIER_IMAG: + break; + case PV_FOURIER_REAL_AND_IMAG: + break; + case PV_FOURIER_PWR: + break; + case PV_FOURIER_PHASE: + break; + default: + break; + } + } else { // not a difference view plot + switch (fCurrentPlotView) { + case PV_DATA: + // get current x-range + xminBin = fMultiGraphData->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fMultiGraphData->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fMultiGraphData->GetXaxis()->GetBinCenter(xminBin); + xmax = fMultiGraphData->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + for (UInt_t i=0; iGetN(); j++) { + // get x and y value + fNonMusrData[i].data->GetPoint(j,xval,yval); + // check if time is in the current range + if ((xval >= xmin) && (xval <= xmax)) { + dump.dataX.push_back(xval); + dump.data.push_back(yval); + dump.dataErr.push_back(fNonMusrData[i].data->GetErrorY(j)); + } + } + + // go through all theory bins + for (Int_t j=0; jGetN(); j++) { + // get x and y value + fNonMusrData[i].theory->GetPoint(j,xval,yval); + // check if time is in the current range + if ((xval >= xmin) && (xval <= xmax)) { + dump.theoryX.push_back(xval); + dump.theory.push_back(yval); + } + } + + // if anything found keep it + if (dump.dataX.size() > 0) + dumpVector.push_back(dump); + } + + break; + case PV_FOURIER_REAL: + break; + case PV_FOURIER_IMAG: + break; + case PV_FOURIER_REAL_AND_IMAG: + break; + case PV_FOURIER_PWR: + break; + case PV_FOURIER_PHASE: + break; + default: + break; + } + } + break; + default: + break; + } + + // generate output filename + + // in order to handle names with "." correctly this slightly odd data-filename generation + TObjArray *tokens = fMsrHandler->GetFileName().Tokenize("."); + TObjString *ostr; + TString str; + TString fln = TString(""); + for (Int_t i=0; iGetEntries()-1; i++) { + ostr = dynamic_cast(tokens->At(i)); + fln += ostr->GetString() + TString("."); + } + if (!fDifferenceView) { + fln += "data.ascii"; + } else { + fln += "diff.ascii"; + } + + if (tokens) { + delete tokens; + tokens = 0; + } + + // open file + ofstream fout; + + // open output data-file + fout.open(fln.Data(), iostream::out); + if (!fout.is_open()) { + cerr << endl << ">> PMusrCanvas::SaveDataAscii: **ERROR** couldn't open file " << fln.Data() << " for writing." << endl; + return; + } + + // find out what is the longest data/theory vector + UInt_t maxDataLength = 0; + UInt_t maxTheoryLength = 0; + for (UInt_t i=0; i 0) + fout << dumpVector[j].dataErr[i] << ", "; + } else { + if (dumpVector[j].dataErr.size() > 0) + fout << ", , , "; + else + fout << ", , "; + } + } + // write last difference entry + if (i 0) + fout << dumpVector[dumpVector.size()-1].dataErr[i]; + } else { + if (dumpVector[dumpVector.size()-1].dataErr.size() > 0) + fout << ", , "; + else + fout << ", "; + } + fout << endl; + } + } else { // no difference view + // write header + switch (fCurrentPlotView) { + case PV_DATA: + fout << "% "; + for (UInt_t i=0; i maxTheoryLength) + maxLength = maxDataLength; + else + maxLength = maxTheoryLength; + + // write data and theory + for (UInt_t i=0; i 0) + fout << dumpVector[j].dataErr[i] << ", "; + } else { + if (dumpVector[j].dataErr.size() > 0) + fout << " , , , "; + else + fout << " , , "; + } + } + // write theory + for (UInt_t j=0; j> Data windows saved in ascii format ..." << endl; + // if (asciiOutput) { + // if (fPlotNumber == static_cast(fMsrHandler->GetMsrPlotList()->size()) - 1) + // Done(0); + // } +} + //-------------------------------------------------------------------------- // CreateStyle (private) //-------------------------------------------------------------------------- @@ -2079,7 +2995,7 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) (Int_t)((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] - data->GetTheoryTimeStart())/data->GetTheoryTimeStep()) * data->GetTheoryTimeStep() - data->GetTheoryTimeStep()/2.0; // closesd start value compatible with the user given end = start + size * data->GetTheoryTimeStep(); // closesd end value compatible with the user given - } +} // invoke histo theoHisto = new TH1F(name, name, size, start, end); @@ -2881,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()); } @@ -2900,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()); } @@ -2919,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()); } @@ -2938,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()); } @@ -2957,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()); } @@ -2976,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()); } @@ -2994,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()); } @@ -3008,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()); } @@ -3022,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()); } @@ -3036,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()); } @@ -3047,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()); } @@ -3069,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()); } @@ -3083,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()); } @@ -3097,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()); } @@ -3111,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()); } @@ -4156,7 +5037,6 @@ void PMusrCanvas::PlotFourier(Bool_t unzoom) // plot fourier data Double_t xmin, xmax, ymin, ymax, binContent; UInt_t noOfPoints = 1000; - switch (fCurrentPlotView) { case PV_FOURIER_REAL: // set x-range @@ -5273,917 +6153,6 @@ void PMusrCanvas::DecrementFourierPhase() } } -//-------------------------------------------------------------------------- -// SaveDataAscii (private) -//-------------------------------------------------------------------------- -/** - *

Saves the currently seen data (data, difference, Fourier spectra, ...) in ascii column format. - */ -void PMusrCanvas::SaveDataAscii() -{ - // collect relevant data - PMusrCanvasAsciiDump dump; - PMusrCanvasAsciiDumpVector dumpVector; - - Int_t xminBin; - Int_t xmaxBin; - Double_t xmin; - Double_t xmax; - Double_t time, freq; - Double_t xval, yval; - - switch (fPlotType) { - case MSR_PLOT_SINGLE_HISTO: - case MSR_PLOT_ASYM: - case MSR_PLOT_MU_MINUS: - if (fDifferenceView) { // difference view plot - switch (fCurrentPlotView) { - case PV_DATA: - // get current x-range - xminBin = fHistoFrame->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fHistoFrame->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fHistoFrame->GetXaxis()->GetBinCenter(xminBin); - xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get time - time = fData[i].diff->GetBinCenter(j); - // check if time is in the current range - if ((time >= xmin) && (time <= xmax)) { - dump.dataX.push_back(time); - dump.data.push_back(fData[i].diff->GetBinContent(j)); - dump.dataErr.push_back(fData[i].diff->GetBinError(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_REAL: - // get current x-range - xminBin = fData[0].diffFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].diffFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_IMAG: - // get current x-range - xminBin = fData[0].diffFourierIm->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].diffFourierIm->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].diffFourierIm->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_REAL_AND_IMAG: - // get current x-range - xminBin = fData[0].diffFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].diffFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].diffFourierRe->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierRe->GetBinContent(j)); - } - } - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_PWR: - // get current x-range - xminBin = fData[0].diffFourierPwr->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].diffFourierPwr->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierPwr->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierPwr->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_PHASE: - // get current x-range - xminBin = fData[0].diffFourierPhase->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].diffFourierPhase->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].diffFourierPhase->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].diffFourierPhase->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - default: - break; - } - } else { // not a difference view plot - switch (fCurrentPlotView) { - case PV_DATA: - // get current x-range - xminBin = fHistoFrame->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fHistoFrame->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fHistoFrame->GetXaxis()->GetBinCenter(xminBin); - xmax = fHistoFrame->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get time - time = fData[i].data->GetBinCenter(j); - // check if time is in the current range - if ((time >= xmin) && (time <= xmax)) { - dump.dataX.push_back(time); - dump.data.push_back(fData[i].data->GetBinContent(j)); - dump.dataErr.push_back(fData[i].data->GetBinError(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get time - time = fData[i].theory->GetBinCenter(j); - // check if time is in the current range - if ((time >= xmin) && (time <= xmax)) { - dump.theoryX.push_back(time); - dump.theory.push_back(fData[i].theory->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - - break; - case PV_FOURIER_REAL: - // get current x-range - xminBin = fData[0].dataFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_IMAG: - // get current x-range - xminBin = fData[0].dataFourierIm->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierIm->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierIm->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_REAL_AND_IMAG: - // get current x-range - xminBin = fData[0].dataFourierRe->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierRe->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierRe->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierRe->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierRe->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - - //----------------------------- - // Im - //----------------------------- - // clean up dump - dump.dataX.clear(); - dump.data.clear(); - dump.dataErr.clear(); - dump.theoryX.clear(); - dump.theory.clear(); - - // go through all data bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierIm->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierIm->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierIm->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - - } - break; - case PV_FOURIER_PWR: - // get current x-range - xminBin = fData[0].dataFourierPwr->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierPwr->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierPwr->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierPwr->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierPwr->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierPwr->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - case PV_FOURIER_PHASE: - // get current x-range - xminBin = fData[0].dataFourierPhase->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierPhase->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetNbinsX(); j++) { - // get frequency - freq = fData[i].dataFourierPhase->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierPhase->GetBinContent(j)); - } - } - - // go through all theory bins - for (Int_t j=1; jGetNbinsX(); j++) { - // get frequency - freq = fData[i].theoryFourierPhase->GetBinCenter(j); - // check if time is in the current range - if ((freq >= xmin) && (freq <= xmax)) { - dump.theoryX.push_back(freq); - dump.theory.push_back(fData[i].theoryFourierPhase->GetBinContent(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - break; - default: - break; - } - } - break; - case MSR_PLOT_NON_MUSR: - if (fDifferenceView) { // difference view plot - switch (fCurrentPlotView) { - case PV_DATA: - // get current x-range - xminBin = fMultiGraphData->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fMultiGraphData->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fMultiGraphData->GetXaxis()->GetBinCenter(xminBin); - xmax = fMultiGraphData->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetN(); j++) { - // get x and y value - fNonMusrData[i].diff->GetPoint(j,xval,yval); - // check if time is in the current range - if ((xval >= xmin) && (xval <= xmax)) { - dump.dataX.push_back(xval); - dump.data.push_back(yval); - dump.dataErr.push_back(fNonMusrData[i].diff->GetErrorY(j)); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - - break; - case PV_FOURIER_REAL: - break; - case PV_FOURIER_IMAG: - break; - case PV_FOURIER_REAL_AND_IMAG: - break; - case PV_FOURIER_PWR: - break; - case PV_FOURIER_PHASE: - break; - default: - break; - } - } else { // not a difference view plot - switch (fCurrentPlotView) { - case PV_DATA: - // get current x-range - xminBin = fMultiGraphData->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fMultiGraphData->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fMultiGraphData->GetXaxis()->GetBinCenter(xminBin); - xmax = fMultiGraphData->GetXaxis()->GetBinCenter(xmaxBin); - - // fill ascii dump data - for (UInt_t i=0; iGetN(); j++) { - // get x and y value - fNonMusrData[i].data->GetPoint(j,xval,yval); - // check if time is in the current range - if ((xval >= xmin) && (xval <= xmax)) { - dump.dataX.push_back(xval); - dump.data.push_back(yval); - dump.dataErr.push_back(fNonMusrData[i].data->GetErrorY(j)); - } - } - - // go through all theory bins - for (Int_t j=0; jGetN(); j++) { - // get x and y value - fNonMusrData[i].theory->GetPoint(j,xval,yval); - // check if time is in the current range - if ((xval >= xmin) && (xval <= xmax)) { - dump.theoryX.push_back(xval); - dump.theory.push_back(yval); - } - } - - // if anything found keep it - if (dump.dataX.size() > 0) - dumpVector.push_back(dump); - } - - break; - case PV_FOURIER_REAL: - break; - case PV_FOURIER_IMAG: - break; - case PV_FOURIER_REAL_AND_IMAG: - break; - case PV_FOURIER_PWR: - break; - case PV_FOURIER_PHASE: - break; - default: - break; - } - } - break; - default: - break; - } - - // generate output filename - - // in order to handle names with "." correctly this slightly odd data-filename generation - TObjArray *tokens = fMsrHandler->GetFileName().Tokenize("."); - TObjString *ostr; - TString str; - TString fln = TString(""); - for (Int_t i=0; iGetEntries()-1; i++) { - ostr = dynamic_cast(tokens->At(i)); - fln += ostr->GetString() + TString("."); - } - if (!fDifferenceView) { - fln += "data.ascii"; - } else { - fln += "diff.ascii"; - } - - if (tokens) { - delete tokens; - tokens = 0; - } - - // open file - ofstream fout; - - // open output data-file - fout.open(fln.Data(), iostream::out); - if (!fout.is_open()) { - cerr << endl << ">> PMusrCanvas::SaveDataAscii: **ERROR** couldn't open file " << fln.Data() << " for writing." << endl; - return; - } - - // find out what is the longest data/theory vector - UInt_t maxDataLength = 0; - UInt_t maxTheoryLength = 0; - for (UInt_t i=0; i 0) - fout << dumpVector[j].dataErr[i] << ", "; - } else { - if (dumpVector[j].dataErr.size() > 0) - fout << ", , , "; - else - fout << ", , "; - } - } - // write last difference entry - if (i 0) - fout << dumpVector[dumpVector.size()-1].dataErr[i]; - } else { - if (dumpVector[dumpVector.size()-1].dataErr.size() > 0) - fout << ", , "; - else - fout << ", "; - } - fout << endl; - } - } else { // no difference view - // write header - switch (fCurrentPlotView) { - case PV_DATA: - fout << "% "; - for (UInt_t i=0; i maxTheoryLength) - maxLength = maxDataLength; - else - maxLength = maxTheoryLength; - - // write data and theory - for (UInt_t i=0; i 0) - fout << dumpVector[j].dataErr[i] << ", "; - } else { - if (dumpVector[j].dataErr.size() > 0) - fout << " , , , "; - else - fout << " , , "; - } - } - // write theory - for (UInt_t j=0; j> Data windows saved in ascii format ..." << endl; -} //-------------------------------------------------------------------------- // IsScaleN0AndBkg (private) @@ -6310,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 new file mode 100644 index 00000000..cf077229 --- /dev/null +++ b/src/classes/PPrepFourier.cpp @@ -0,0 +1,438 @@ +/*************************************************************************** + + PPrepFourier.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 "PPrepFourier.h" + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor. + */ +PPrepFourier::PPrepFourier() +{ + fBkgRange[0] = -1; + fBkgRange[1] = -1; + fPacking = 1; +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor. + */ +PPrepFourier::PPrepFourier(const Int_t *bkgRange, const Int_t packing) : + fPacking(packing) +{ + SetBkgRange(bkgRange); +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor. + */ +PPrepFourier::~PPrepFourier() +{ + fRawData.clear(); + fData.clear(); +} + +//-------------------------------------------------------------------------- +// SetBkgRange +//-------------------------------------------------------------------------- +/** + *

set the background range. + * + * \param bkgRange array with background range + */ +void PPrepFourier::SetBkgRange(const Int_t *bkgRange) +{ + int err=0; + if (bkgRange[0] >= -1) { + fBkgRange[0] = bkgRange[0]; + } else { + err = 1; + } + if (bkgRange[1] >= -1) { + fBkgRange[1] = bkgRange[1]; + } else { + if (err == 0) + err = 2; + else + err = 3; + } + + TString errMsg(""); + switch (err) { + case 1: + errMsg = TString::Format("start bkg range < 0 (given: %d), will ignore it.", bkgRange[0]); + break; + case 2: + errMsg = TString::Format("end bkg range < 0 (given: %d), will ignore it.", bkgRange[1]); + break; + case 3: + errMsg = TString::Format("start/end bkg range < 0 (given: %d/%d), will ignore it.", bkgRange[0], bkgRange[1]); + break; + default: + errMsg = TString("??"); + break; + } + + if (err != 0) { + cerr << endl << ">> PPrepFourier::SetBkgRange: **WARNING** " << errMsg << endl; + } +} + +//-------------------------------------------------------------------------- +// SetPacking +//-------------------------------------------------------------------------- +/** + *

set the packing for the histograms. + * + * \param packing number to be used. + */ +void PPrepFourier::SetPacking(const Int_t packing) +{ + if (packing > 0) { + fPacking = packing; + } else { + cerr << endl << ">> PPrepFourier::SetPacking: **WARNING** found packing=" << packing << " < 0, will ignore it." << endl; + } +} + +//-------------------------------------------------------------------------- +// AddData +//-------------------------------------------------------------------------- +/** + *

add a data-set (time domain data + meta information) to the internal + * data vector. + * + * \param data set to be added + */ +void PPrepFourier::AddData(musrFT_data &data) +{ + fRawData.push_back(data); +} + +//-------------------------------------------------------------------------- +// DoBkgCorrection +//-------------------------------------------------------------------------- +/** + *

Correct the internal data sets according to a background interval given. + */ +void PPrepFourier::DoBkgCorrection() +{ + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } + + // if no bkg-range is given, nothing needs to be done + if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) { + return; + } + + // make sure that the bkg range is ok + for (unsigned int i=0; i= fRawData[i].rawData.size()) || (fBkgRange[1] >= fRawData[i].rawData.size())) { + cerr << endl << "PPrepFourier::DoBkgCorrection() **ERROR** bkg-range out of data-range!"; + return; + } + } + + Double_t bkg=0.0; + for (unsigned int i=0; iRebin (pack) the internal data. + */ +void PPrepFourier::DoPacking() +{ + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } + + if (fPacking == 1) // nothing to be done + return; + + PDoubleVector tmpData; + Double_t dval = 0.0; + for (unsigned int i=0; iNot implemented yet. + */ +void PPrepFourier::DoFiltering() +{ + // 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 rescaling factor for the estimated N0. Should be around 1 + */ +void PPrepFourier::DoLifeTimeCorrection(Double_t fudge) +{ + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } + + // calc exp(+t/tau)*N(t), where N(t) is already background corrected + Double_t scale; + for (unsigned int i=0; iReturns the meta information of a data set. + * + * \param idx index of the object + */ +TString PPrepFourier::GetInfo(const UInt_t idx) +{ + TString info(""); + + if (idx < fRawData.size()) + info = fRawData[idx].info; + + return info; +} + +//-------------------------------------------------------------------------- +// GetData +//-------------------------------------------------------------------------- +/** + *

Creates the requested TH1F objects and returns them. The ownership is with + * the caller. + */ +vector PPrepFourier::GetData() +{ + vector data; + data.resize(fData.size()); + + // if not data are present, just return an empty vector + if (fData.size() == 0) + return data; + + TString name(""); + Double_t dt=0.0; + Double_t start=0.0; + Double_t end=0.0; + UInt_t size; + UInt_t startIdx; + UInt_t endIdx; + + for (unsigned int i=0; i= 0.0) { + startIdx = (UInt_t)(start/dt)+1; + endIdx = (UInt_t)(end/dt)+1; + } else { + cerr << endl << ">> PPrepFourier::GetData **WARNING** found start time < 0.0, will set it to 0.0" << endl; + endIdx = static_cast(end/dt)+1; + } + } + + if (start == -1.0) { // no time range given, hence start == 0.0 - dt/2 + start = -dt/2.0; + } else { // time range given + start -= dt/2.0; + } + + if (end == -1.0) { // no time range given, hence end == (fData[idx].size()-1)*dt + dt/2 + end = (fData[i].size()-1)*dt+dt/2.0; + } else { // time range given + end += dt/2.0; + } + + data[i] = new TH1F(name.Data(), fRawData[i].info.Data(), size, start, end); + for (unsigned int j=startIdx; jSetBinContent(j-startIdx+1, fData[i][j]); + } + + return data; +} + +//-------------------------------------------------------------------------- +// GetData +//-------------------------------------------------------------------------- +/** + *

Creates the requested TH1F object and returns it. The ownership is with + * the caller. + * + * \param idx index of the requested histogram + */ +TH1F *PPrepFourier::GetData(const UInt_t idx) +{ + if (fData.size() == 0) // no data present + return 0; + + if (idx > fData.size()) // requested index out of range + return 0; + + TString name = TString::Format("histo%2d", idx); + Double_t dt = fRawData[idx].timeResolution*fPacking; + Double_t start = fRawData[idx].timeRange[0]; + Double_t end = fRawData[idx].timeRange[1]; + UInt_t size = fData[idx].size(); + UInt_t startIdx = 1; + UInt_t endIdx = size; + + // time range given, hence calculate the proper size + if (start != -1.0) { + size = (UInt_t)((end-start)/dt); + if (start >= 0.0) { + startIdx = (UInt_t)(start/dt)+1; + endIdx = (UInt_t)(end/dt)+1; + } else { + cerr << endl << ">> PPrepFourier::GetData **WARNING** found start time < 0.0, will set it to 0.0" << endl; + endIdx = (UInt_t)(end/dt)+1; + } + } + + if (start == -1.0) { // no time range given, hence start == 0.0 - dt/2 + start = -dt/2.0; + } else { // time range given + start -= dt/2.0; + } + + if (end == -1.0) { // no time range given, hence end == (fData[idx].size()-1)*dt + dt/2 + end = (fData[idx].size()-1)*dt+dt/2.0; + } else { // time range given + end += dt/2.0; + } + + TH1F *data = new TH1F(name.Data(), fRawData[idx].info.Data(), size, start, end); + for (unsigned int i=startIdx; iSetBinContent(i-startIdx+1, fData[idx][i]); + + return data; +} + +//-------------------------------------------------------------------------- +// InitData +//-------------------------------------------------------------------------- +/** + *

Copy raw-data to internal data from t0 to the size of raw-data. + */ +void PPrepFourier::InitData() +{ + fData.resize(fRawData.size()); + unsigned int t0; + for (unsigned int i=0; i= 0) + t0 = fRawData[i].t0; + else + t0 = 0; + for (unsigned int j=t0; jChecks if runName is found, and if so return these data. + * runName is as given in the msr-file. * * return: * - if data are found: pointer to the data. @@ -240,6 +241,26 @@ PRawRunData* PRunDataHandler::GetRunData(const TString &runName) return &fData[i]; } +//-------------------------------------------------------------------------- +// GetRunData +//-------------------------------------------------------------------------- +/** + *

return data-set with index idx. + * + * return: + * - if data are found: pointer to the data. + * - otherwise the null pointer will be returned. + * + * \param idx index of the raw data set. + */ +PRawRunData* PRunDataHandler::GetRunData(const UInt_t idx) +{ + if (idx >= fData.size()) + return 0; + else + return &fData[idx]; +} + //-------------------------------------------------------------------------- // ReadData //-------------------------------------------------------------------------- diff --git a/src/external/libBNMR/045674.msr b/src/external/libBNMR/045674.msr new file mode 100644 index 00000000..df6bc823 Binary files /dev/null and b/src/external/libBNMR/045674.msr differ diff --git a/src/external/libBNMR/ExpRlx-MUD.msr b/src/external/libBNMR/ExpRlx-MUD.msr new file mode 100644 index 00000000..ccbf0bb4 --- /dev/null +++ b/src/external/libBNMR/ExpRlx-MUD.msr @@ -0,0 +1,73 @@ + +# Run Numbers: 1111 +############################################################### +FITPARAMETER +############################################################### +# No Name Value Err Min Max + 1 Alphap 1.11662 0.00020 none + 2 Asyp 0 0.00038 none + 3 T 1e6 0 none + 4 Rlx 0.969e6 0.020 none + 5 Pos 1 0 none + 6 Neg -1 0 none + +############################################################### +THEORY +############################################################### +asymmetry fun1 +userFcn libBNMR.so ExpRlx 3 4 + +############################################################### +FUNCTIONS +############################################################### +fun1 = map1 * map2 + +############################################################### +RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +forward 3 +backward 4 +data 11 799 11 799 +background 800 900 800 900 # estimated bkg: 416.9700 / 465.7600 +t0 0.0 0.0 +map 2 5 0 0 0 0 0 0 0 0 +fit 5e5 8e6 +packing 5 + +RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +forward 5 +backward 6 +data 11 799 11 799 +background 800 900 800 900 # estimated bkg: 430.9200 / 479.4500 +t0 0.0 0.0 +map 2 6 0 0 0 0 0 0 0 0 +fit 5e5 8e6 +packing 5 + + +############################################################### +COMMANDS +MINIMIZE +HESSE +SAVE + +############################################################### +PLOT 2 (asymmetry plot) +runs 1 2 +use_fit_ranges + + +############################################################### +FOURIER +units MHz # units either 'Gauss', 'MHz', or 'Mc/s' +fourier_power 12 +apodization STRONG # NONE, WEAK, MEDIUM, STRONG +plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE +phase 8 +#range FRQMIN FRQMAX +############################################################### +STATISTIC --- 2014-11-14 16:37:41 + chisq = 372.2, NDF = 295, chisq/NDF = 1.261744 diff --git a/src/include/PFourier.h b/src/include/PFourier.h index e8355290..787450e1 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -8,7 +8,7 @@ ***************************************************************************/ /*************************************************************************** - * Copyright (C) 2007-2014 by Andreas Suter * + * Copyright (C) 2007-2015 by Andreas Suter * * andreas.suter@psi.ch * * * * This program is free software; you can redistribute it and/or modify * @@ -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); @@ -68,7 +70,7 @@ class PFourier Int_t fApodization; ///< 0=none, 1=weak, 2=medium, 3=strong - Double_t fTimeResolution; ///< time resolution of the data histogram + Double_t fTimeResolution; ///< time resolution of the data histogram in (us) Double_t fStartTime; ///< start time of the data histogram Double_t fEndTime; ///< end time of the data histogram Bool_t fDCCorrected; ///< if true, removed DC offset from signal before Fourier transformation, otherwise not 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/PMusr.h b/src/include/PMusr.h index f9826fbd..f083b13b 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -8,7 +8,7 @@ ***************************************************************************/ /*************************************************************************** - * Copyright (C) 2007-2014 by Andreas Suter * + * Copyright (C) 2007-2015 by Andreas Suter * * andreas.suter@psi.ch * * * * This program is free software; you can redistribute it and/or modify * @@ -44,6 +44,7 @@ using namespace std; #define PMUSR_TOKENIZE_ERROR -5 #define PMUSR_MSR_LOG_FILE_WRITE_ERROR -6 #define PMUSR_MSR_FILE_WRITE_ERROR -7 +#define PMUSR_DATA_FILE_READ_ERROR -8 #define PRUN_SINGLE_HISTO 0 #define PRUN_ASYMMETRY 2 @@ -703,7 +704,7 @@ typedef vector PMsrRunList; */ typedef struct { Bool_t fFourierBlockPresent; ///< flag indicating if a Fourier block is present in the msr-file - Int_t fUnits; ///< flag used to indicate the units. 0=field units (G); 1=frequency units (MHz); 2=Mc/s + Int_t fUnits; ///< flag used to indicate the units. 1=field units (G); 2=field units (T); 3=frequency units (MHz); 4=Mc/s Bool_t fDCCorrected; ///< if set true, the dc offset of the signal/theory will be removed before the FFT is made. Int_t fFourierPower; ///< i.e. zero padding up to 2^fFourierPower, default = 0 which means NO zero padding Int_t fApodization; ///< tag indicating the kind of apodization wished, 0=no appodization (default), 1=weak, 2=medium, 3=strong (for details see the docu) diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index eaf0b5ac..e8296e4c 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -229,6 +229,7 @@ class PMusrCanvas : public TObject, public TQObject virtual void LastCanvasClosed(); // SLOT virtual void SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat); + virtual void SaveDataAscii(); private: Int_t fTimeout; ///< timeout after which the Done signal should be emited. If timeout <= 0, no timeout is taking place @@ -277,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) @@ -329,12 +328,10 @@ class PMusrCanvas : public TObject, public TQObject virtual void IncrementFourierPhase(); virtual void DecrementFourierPhase(); - virtual void SaveDataAscii(); - 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 new file mode 100644 index 00000000..d1066809 --- /dev/null +++ b/src/include/PPrepFourier.h @@ -0,0 +1,88 @@ +/*************************************************************************** + + PPrepFourier.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. * + ***************************************************************************/ + +#ifndef _PPREPFOURIER_H_ +#define _PPREPFOURIER_H_ + +#include +#include +using namespace std; + +#include "TH1F.h" + +#include "PMusr.h" + +//---------------------------------------------------------------------------- +/** + *

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 +} musrFT_data; + +//---------------------------------------------------------------------------- +/** + *

Little helper class to prepare time-domain data for Fourier transform, without + * theory, etc. + */ +class PPrepFourier { + public: + PPrepFourier(); + PPrepFourier(const Int_t *bkgRange, const Int_t packing); + virtual ~PPrepFourier(); + + void SetBkgRange(const Int_t *bkgRange); + void SetPacking(const Int_t packing); + void AddData(musrFT_data &data); + void DoBkgCorrection(); + void DoPacking(); + void DoFiltering(); + void DoLifeTimeCorrection(Double_t fudge); + + TString GetInfo(const UInt_t idx); + UInt_t GetNoOfData() { return fRawData.size(); } + vector GetData(); + TH1F *GetData(const UInt_t idx); + + private: + vector fRawData; + vectorfData; + Int_t fBkgRange[2]; + Int_t fPacking; + + void InitData(); +}; + +#endif // _PPREPFOURIER_H_ + diff --git a/src/include/PRunDataHandler.h b/src/include/PRunDataHandler.h index ca7e4e79..68c27eab 100644 --- a/src/include/PRunDataHandler.h +++ b/src/include/PRunDataHandler.h @@ -8,7 +8,7 @@ ***************************************************************************/ /*************************************************************************** - * Copyright (C) 2007-2014 by Andreas Suter * + * Copyright (C) 2007-2015 by Andreas Suter * * andreas.suter@psi.ch * * * * This program is free software; you can redistribute it and/or modify * @@ -60,6 +60,7 @@ class PRunDataHandler virtual Bool_t IsAllDataAvailable() const { return fAllDataAvailable; } virtual PRawRunData* GetRunData(const TString &runName); + virtual PRawRunData* GetRunData(const UInt_t idx=0); virtual TString GetRunPathName() {return fRunPathName; } private: diff --git a/src/musrFT.cpp b/src/musrFT.cpp new file mode 100644 index 00000000..68d3c47b --- /dev/null +++ b/src/musrFT.cpp @@ -0,0 +1,1428 @@ +/*************************************************************************** + + musrFT.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. * + ***************************************************************************/ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include + +#include +#include +#include +using namespace std; + +#include +#include +#include +#include +#include +#include + +#include "git-revision.h" +#include "PMusr.h" +#include "PStartupHandler.h" +#include "PMsrHandler.h" +#include "PRunDataHandler.h" +#include "PPrepFourier.h" +#include "PFourier.h" +#include "PFourierCanvas.h" + +//---------------------------------------------------------------------------- +/** + *

Structure keeping the command line options. + */ +typedef struct { + 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() +{ + cout << endl << "usage: musrFT [Options] [ | -df, --data-file ]"; + cout << endl << " : msr-file name(s). These msr-files are used for the Fourier transform."; + cout << endl << " It can be a list of msr-files, e.g. musrFT 3110.msr 3111.msr"; + cout << endl << " For the syntax of the msr-file check the user manual of musrfit."; + cout << endl << " -df, --data-file : This allows to feed only muSR data file(s) to"; + cout << endl << " perform the Fourier transform. Since the extended information"; + cout << endl << " are missing, they will need to be provided by to options, or musrFT"; + cout << endl << " tries to guess, based on musrFT_startup.xml settings."; + cout << endl << " Options: "; + cout << endl << " --help : display this help and exit"; + cout << endl << " --version : output version information and exit"; + cout << endl << " -g, --graphic-format : "; + cout << endl << " will produce a graphic-output-file without starting a root session."; + cout << endl << " the name is based either on the or the ,"; + cout << endl << " e.g. 3310.msr -> 3310_0.png."; + cout << endl << " Supported graphic-format-extension: eps, pdf, gif, jpg, png, svg, xpm, root"; + cout << endl << " --dump : rather than starting a root session and showing Fourier graphs of the data,"; + cout << endl << " it will output the Fourier data in an ascii file ."; + 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', '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 '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."; + cout << endl << " Default is -1, i.e. no zero padding will be performed."; + cout << endl << " -u, --units : is used to define the x-axis of the Fourier transform."; + 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 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."; + cout << endl << " Default will be -1.0 for both which means, take the full Fourier range."; + cout << endl << " -tr, --time-range : time domain range to be used for Fourier transform."; + cout << endl << " , are to be given in (us). If nothing is given, the full time range"; + cout << endl << " found in the data file(s) will be used."; + cout << endl << " --histo : give the of histograms to be used for the Fourier transform."; + cout << endl << " E.g. musrFT -df lem15_his_01234.root --histo 1 3, will only be needed together with"; + cout << endl << " the option --data-file. If multiple data file are given, will apply"; + cout << endl << " to all data-files given. If --histo is not given, all histos of a data file will be used."; + cout << endl << " -a, --average : show the average of all Fourier transformed data."; + cout << endl << " --t0 : A list of t0's can be provided. This in conjunction with --data-file and"; + cout << endl << " --fourier-option real allows to get the proper inital phase if t0's are known."; + cout << endl << " If a single t0 for multiple histos is given, it is assume, that this t0 is common"; + 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 << " --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 << " -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; +} + +//------------------------------------------------------------------------- +/** + * <p>initialize startup parameters. + * + * \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("??"); + startupParam.apodization = TString("none"); + startupParam.fourierPower = -1; + startupParam.fourierUnits = TString("??"); + startupParam.initialPhase = 0.0; + startupParam.fourierRange[0] = -1.0; + startupParam.fourierRange[1] = -1.0; + startupParam.timeRange[0] = -1.0; + startupParam.timeRange[1] = -1.0; + startupParam.showAverage = false; + startupParam.packing = 1; + startupParam.title = TString(""); + startupParam.lifetimecorrection = 0.0; + startupParam.timeout = 3600; +} + +//------------------------------------------------------------------------- +/** + * <p>Parses the musrFT --histo options. Allowed --histo options are: + * \<h0\> \<h1\> ... \<hN\>, e.g. --histo 1 3 5 + * \<h0\>-\<hN\>, e.g. --histo 1-16 + * or a combination of both + * e.g. --histo 1 3 7-12 15 + * + * <b>return:</b> true if everything is OK, false otherwise. + * + * \param i position of the --histo option within argv. At return it will be shifted + * to the last element of the --histo option. + * \param argc number of elements in argv + * \param argv list of command line tokens + * \param startupParam startup parameter structure + */ +bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &startupParam) +{ + int start = i+1, end = 0; + + // find last element of histo option + while (++i < argc) { + if (argv[i][0] == '-') { + if (!isdigit(argv[i][1]) || (argv[i][1] == '-')) + break; + } + } + end = i; + --i; + if (end < start) { + cerr << endl << ">> musrFT **ERROR** something is wrong with the --histo arguments." << endl; + startupParam.histo.clear(); + return false; + } + + // handle histo arguments + TString tstr(""); + for (int j=start; j<end; j++) { + tstr = argv[j]; + if (!tstr.Contains("-")) { // likely to be a single number + if (tstr.IsDigit()) { + startupParam.histo.push_back(tstr.Atoi()); + } else { // not a number -> error + cerr << endl << ">> musrFT **ERROR** found --histo argument '" << tstr << "' which is not a number." << endl; + startupParam.histo.clear(); + return false; + } + } else { // should be something like h0-hN with h0, hN numbers + TObjArray *tok = tstr.Tokenize("-"); + if (tok->GetEntries() != 2) { + cerr << endl << ">> musrFT **ERROR** found --histo argument '" << tstr << "' which is not of the form <h0>-<hN>." << endl; + startupParam.histo.clear(); + return false; + } + TObjString *ostr; + TString sstr(""); + int first=0, last=0; + ostr = dynamic_cast<TObjString*>(tok->At(0)); + sstr = ostr->GetString(); + if (sstr.IsDigit()) { + first = sstr.Atoi(); + } else { + cerr << endl << ">> musrFT **ERROR** found --histo argument '" << tstr << "' which is of the form <h0>-<hN>,"; + cerr << endl << " but <h0>='" << sstr << "' is not a number." << endl; + startupParam.histo.clear(); + return false; + } + ostr = dynamic_cast<TObjString*>(tok->At(1)); + sstr = ostr->GetString(); + if (sstr.IsDigit()) { + last = sstr.Atoi(); + } else { + cerr << endl << ">> musrFT **ERROR** found --histo argument '" << tstr << "' which is of the form <h0>-<hN>,"; + cerr << endl << " but <hN>='" << sstr << "' is not a number." << endl; + startupParam.histo.clear(); + return false; + } + + if (first > last) { + cerr << endl << ">> musrFT **ERROR** found --histo argument of the form <h0>-<hN> with h0=" << first << " > hN=" << last << "." << endl; + startupParam.histo.clear(); + return false; + } + + for (int k=first; k<=last; k++) { + startupParam.histo.push_back(k); + } + + // clean up + if (tok) + delete tok; + } + } + + return true; +} + +//------------------------------------------------------------------------- +/** + * <p>Parses the musrFT command line options. + * + * <b>return:</b> 0 if everything is OK, 1 for --version or --help, 2 for an error. + * + * \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) +{ + TString tstr(""); + + for (int i=1; i<argc; i++) { + tstr = argv[i]; + if (tstr.BeginsWith("--version")) { +#ifdef HAVE_CONFIG_H + cout << endl << "musrFT version: " << PACKAGE_VERSION << ", git-rev: " << GIT_REVISION << endl << endl; +#else + cout << endl << "musrFT git-rev: " << GIT_REVISION << endl << endl; +#endif + return 1; + } else if (tstr.BeginsWith("--help")) { + musrFT_syntax(); + return 1; + } else if (tstr.BeginsWith("-g") || tstr.BeginsWith("--graphic-format")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --graphic-format without argument!" << endl; + return 2; + } + TString topt(argv[i+1]); + if (!topt.BeginsWith("eps") && !topt.BeginsWith("pdf") && !topt.BeginsWith("gif") && !topt.BeginsWith("jpg") && + !topt.BeginsWith("png") && !topt.BeginsWith("svg") && !topt.BeginsWith("xpm") && !topt.BeginsWith("root")) { + cerr << endl << ">> musrFT **ERROR** found unrecogniced graphic format '" << topt << "'!" << endl; + return 2; + } + startupParam.graphicFormat = topt; + i++; + } else if (tstr.BeginsWith("--dump")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --dump without argument!" << endl; + return 2; + } + startupParam.dumpFln = argv[i+1]; + i++; + } else if (tstr.Contains("--filter")) { + cout << endl << "debug> found option filter. NOT YET ANY FUNCTIONALITY." << endl; + } else if (tstr.Contains("-b") || tstr.Contains("--background")) { + if (i+2 >= argc) { // something is wrong since there needs to be two arguments here + cerr << endl << ">> musrFT **ERROR** found option --background with wrong number of arguments." << endl; + return 2; + } + TString bkg[2]; + bkg[0] = argv[i+1]; + bkg[1] = argv[i+2]; + if (!bkg[0].IsDigit()) { + cerr << endl << ">> musrFT **ERROR** <start> bin of option --background is NOT an int-number! ('" << bkg[0] << "')." << endl; + return 2; + } + if (!bkg[1].IsDigit()) { + cerr << endl << ">> musrFT **ERROR** <end> bin of option --background is NOT an int-number! ('" << bkg[1] << "')." << endl; + return 2; + } + startupParam.bkg[0] = bkg[0].Atoi(); + startupParam.bkg[1] = bkg[1].Atoi(); + i += 2; + } else if (tstr.BeginsWith("-fo") || tstr.BeginsWith("--fourier-option")) { + if (i+1 >= argc) { // something is wrong since there needs to be two arguments here + cerr << endl << ">> musrFT **ERROR** found option --fourier-option without arguments." << endl; + return 2; + } + TString topt(argv[i+1]); + if (!topt.BeginsWith("real") && !topt.BeginsWith("imag") && !topt.BeginsWith("power") && !topt.BeginsWith("phase")) { + cerr << endl << ">> musrFT **ERROR** found option --fourier-option with unrecognized argument '" << topt << "'." << endl; + return 2; + } + startupParam.fourierOpt = topt; + i++; + } else if (tstr.BeginsWith("-apod") || tstr.BeginsWith("--apodization")) { + if (i+1 >= argc) { // something is wrong since there needs to be two arguments here + cerr << endl << ">> musrFT **ERROR** found option --apodization without arguments." << endl; + return 2; + } + TString topt(argv[i+1]); + if (!topt.BeginsWith("none") && !topt.BeginsWith("weak") && !topt.BeginsWith("medium") && !topt.BeginsWith("strong")) { + cerr << endl << ">> musrFT **ERROR** found option --apodization with unrecognized argument '" << topt << "'." << endl; + return 2; + } + 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 + cerr << endl << ">> musrFT **ERROR** found option --fourier-power without arguments." << endl; + return 2; + } + TString fourierPower(argv[i+1]); + if (!fourierPower.IsDigit()) { + cerr << endl << ">> musrFT **ERROR** found option --fourier-power with a power which is not an integer '" << fourierPower << "'." << endl; + return 2; + } + startupParam.fourierPower = fourierPower.Atoi(); + if ((startupParam.fourierPower < 1) || (startupParam.fourierPower > 20)) { + cerr << endl << ">> musrFT **ERROR** found Fourier power '" << fourierPower << "', which is out of range [1..20]" << endl; + return 2; + } + i++; + } else if (tstr.BeginsWith("-u") || tstr.BeginsWith("--units")) { + if (i+1 >= argc) { // something is wrong since there needs to be two arguments here + cerr << endl << ">> musrFT **ERROR** found option --units without arguments." << endl; + return 2; + } + TString topt(argv[i+1]); + if (!topt.BeginsWith("MHz", TString::kIgnoreCase) && !topt.BeginsWith("Gauss", TString::kIgnoreCase) && + !topt.BeginsWith("Tesla", TString::kIgnoreCase) && !topt.BeginsWith("Mc/s", TString::kIgnoreCase)) { + cerr << endl << ">> musrFT **ERROR** found option --fourier-option with unrecognized argument '" << topt << "'." << endl; + return 2; + } + startupParam.fourierUnits = topt; + i++; + } else if (tstr.BeginsWith("-ph") || tstr.BeginsWith("--phase")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --phase without argument!" << endl; + return 2; + } + TString phase(argv[i+1]); + if (!phase.IsFloat()) { + cerr << endl << ">> musrFT **ERROR** found --phase argument '" << phase << "' which is not a number." << endl; + return 2; + } + startupParam.initialPhase = phase.Atof(); + i++; + } else if (tstr.BeginsWith("-fr") || tstr.BeginsWith("--fourier-range")) { + if (i+2 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --fourier-range with wrong number of arguments!" << endl; + return 2; + } + TString fourierRange[2] = {argv[i+1], argv[i+2]}; + if (!fourierRange[0].IsFloat() || !fourierRange[1].IsFloat()) { + cerr << endl << ">> musrFT **ERROR** found invalid --fourier-range arguments '" << fourierRange[0] << "' and/or '" << fourierRange[1] << "'." << endl; + return 2; + } + startupParam.fourierRange[0] = fourierRange[0].Atof(); + startupParam.fourierRange[1] = fourierRange[1].Atof(); + i += 2; + } else if (tstr.BeginsWith("-tr") || tstr.BeginsWith("--time-range")) { + if (i+2 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --time-range with wrong number of arguments!" << endl; + return 2; + } + TString timeRange[2] = {argv[i+1], argv[i+2]}; + if (!timeRange[0].IsFloat() || !timeRange[1].IsFloat()) { + cerr << endl << ">> musrFT **ERROR** found invalid --time-range arguments '" << timeRange[0] << "' and/or '" << timeRange[1] << "'." << endl; + return 2; + } + startupParam.timeRange[0] = timeRange[0].Atof(); + startupParam.timeRange[1] = timeRange[1].Atof(); + i += 2; + } else if (tstr.BeginsWith("-a") || tstr.BeginsWith("--average")) { + startupParam.showAverage = true; + } else if (tstr.BeginsWith("--histo")) { + if (!musrFT_filter_histo(i, argc, argv, startupParam)) + return 2; + } else if (tstr.BeginsWith("--t0")) { + TString topt(""); + while (++i < argc) { + if (argv[i][0] == '-') { + --i; + break; + } else { + topt = argv[i]; + if (!topt.IsDigit()) { + cerr << endl << ">> musrFT **ERROR** found option t0='" << topt << "' which is not a number" << endl; + return 2; + } + startupParam.t0.push_back(topt.Atoi()); + } + } + if (startupParam.dataFln.size() == 0) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --t0 without argument!" << endl; + return 2; + } + } 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; + } + ++i; + startupParam.title = argv[i]; + } else if (tstr.BeginsWith("-pa") || tstr.BeginsWith("--packing")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --packing without argument!" << endl; + return 2; + } + ++i; + TString pack = TString(argv[i]); + if (!pack.IsDigit()) { + cerr << endl << ">> musrFT **ERROR** found option --packing with argument '" << pack << "' which is NOT an integer!" << endl; + 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; + return 2; + } + ++i; + TString fudge(argv[i]); + if (!fudge.IsFloat()) { + 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] == '-') { + --i; + break; + } else { + startupParam.dataFln.push_back(argv[i]); + TString fln(argv[i]); + TString fileFormat("??"); + if (fln.Contains(".root", TString::kIgnoreCase)) + fileFormat = "MusrRoot"; + else if (fln.Contains(".bin", TString::kIgnoreCase) || fln.Contains(".mdu", TString::kIgnoreCase)) + fileFormat = "PsiBin"; + else if (fln.Contains(".nxs", TString::kIgnoreCase)) + fileFormat = "NeXus"; + else if (fln.Contains(".msr", TString::kIgnoreCase)) + fileFormat = "Mud"; + + if (fileFormat == "??") { + cerr << endl << ">> musrFT **ERROR** found data file name with unrecognized data file format ('" << fln << "')." << endl; + return 2; + } else { + startupParam.dataFileFormat.push_back(fileFormat); + } + } + } + if (startupParam.dataFln.size() == 0) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --data-file without argument!" << endl; + return 2; + } + } else if (tstr.Contains(".msr")) { + startupParam.msrFln.push_back(tstr); + } else { + cerr << endl << ">> musrFT **ERROR** unrecognized option '" << tstr << "' found." << endl; + return 2; + } + } + + // consistency checks + if ((startupParam.msrFln.size() == 0) && (startupParam.dataFln.size() == 0)) { + cerr << endl << ">> musrFT **ERROR** neither <msr-file> nor <data-file> defined." << endl; + return 2; + } + if (startupParam.bkg[0] > startupParam.bkg[1]) { + cerr << endl << ">> musrFT **WARNING** in --background, start=" << startupParam.bkg[0] << " > end=" << startupParam.bkg[1] << ", will swap them." << endl; + double swap = startupParam.bkg[0]; + startupParam.bkg[0] = startupParam.bkg[1]; + startupParam.bkg[1] = swap; + } + if (startupParam.fourierRange[0] > startupParam.fourierRange[1]) { + cerr << endl << ">> musrFT **WARNING** in --fourier-range, start=" << startupParam.fourierRange[0] << " > end=" << startupParam.fourierRange[1] << ", will swap them." << endl; + double swap = startupParam.fourierRange[0]; + startupParam.fourierRange[0] = startupParam.fourierRange[1]; + startupParam.fourierRange[1] = swap; + } + + return 0; +} + +//---------------------------------------------------------------------------------------- +/** + * <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) +{ + double dval; + TString str = fln; + // file name + // trunc it in case a path-name is given + size_t idx = str.Last('/'); + if (idx > 0) + str.Remove(0, idx+1); + metaInfo = str; + metaInfo += ","; + // temperature + for (UInt_t i=0; i<rawRunData->GetNoOfTemperatures(); i++) { + metaInfo += TString::Format("T%d=%0.2fK,", i, rawRunData->GetTemperature(i)); + } + // magnetic field + dval = rawRunData->GetField(); + if (dval == PMUSR_UNDEFINED) + metaInfo += TString("B=??,"); + else if (dval < 5000.0) + metaInfo += TString::Format("B=%0.1fG,", dval); + else + metaInfo += TString::Format("B=%0.1fT,", dval/1.0e4); + // implantation energy + dval = rawRunData->GetEnergy(); + if (dval == PMUSR_UNDEFINED) + metaInfo += TString("E=??;"); + else if (dval < 1000.0) + metaInfo += TString::Format("E=%0.1fkeV;", dval); + else + metaInfo += TString::Format("E=%0.1fMeV;", dval/1.0e3); + + metaInfo += *rawRunData->GetCryoName(); + metaInfo += ";"; + metaInfo += *rawRunData->GetSample(); +} + +//------------------------------------------------------------------------- +/** + * <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"; + cout << endl << " '" << rd.info << "'"; + cout << endl << " NO warranty this is sensible!" << endl; + + unsigned int idx = 0; + double max = rd.rawData[0]; + for (unsigned int i=1; i<rd.rawData.size(); i++) { + if (rd.rawData[i] > max) { + max = rd.rawData[i]; + idx = (int)i; + } + } + cout << endl << ">> musrFT_estimateT0: estimated t0=" << idx << endl; + rd.t0 = idx; +} + +//------------------------------------------------------------------------- +/** + * <p> deletes a histogram. + * + * \param h point to a ROOT histogram object + */ +void musrFT_cleanup(TH1F *h) +{ + if (h) { + delete h; + h = 0; + } +} + +//------------------------------------------------------------------------- +/** + * <p>Dump the Fourier transformed data into an ascii file. + * + * \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) +{ + vector<PDoubleVector> data; + PDoubleVector freq; + PDoubleVector re; + PDoubleVector im; + PDoubleVector pwr; + TH1F *hRe=0, *hIm=0; + Double_t dval; + + // make sure start/end are given, otherwise take the minimum/maximum off all data + hRe = fourierData[0]->GetRealFourier(); + if (start == -1.0) { + start = hRe->GetBinCenter(1); + if (end == -1.0) + end = hRe->GetBinCenter(hRe->GetNbinsX()); + } + + unsigned int minSize = hRe->GetNbinsX()-1; + musrFT_cleanup(hRe); + for (unsigned int i=1; i<fourierData.size(); i++) { + hRe = fourierData[i]->GetRealFourier(); + if (hRe->GetNbinsX()-1 < minSize) + minSize = hRe->GetNbinsX()-1; + musrFT_cleanup(hRe); + } + + for (unsigned int i=0; i<fourierData.size(); i++) { + hRe = fourierData[i]->GetRealFourier(); + hIm = fourierData[i]->GetImaginaryFourier(); + for (int j=1; j<minSize; j++) { + dval = hRe->GetBinCenter(j); + if ((dval >= start) && (dval <= end)) { + freq.push_back(dval); + re.push_back(hRe->GetBinContent(j)); + im.push_back(hIm->GetBinContent(j)); + pwr.push_back(hRe->GetBinContent(j)*hRe->GetBinContent(j)+hIm->GetBinContent(j)*hIm->GetBinContent(j)); + } + } + data.push_back(freq); + data.push_back(re); + data.push_back(im); + data.push_back(pwr); + // cleanup + freq.clear(); + re.clear(); + im.clear(); + pwr.clear(); + musrFT_cleanup(hRe); + musrFT_cleanup(hIm); + } + + ofstream fout(fln, ofstream::out); + + // write header + fout << "% "; + for (unsigned int i=0; i<fourierData.size()-1; i++) + fout << "freq" << i << ", Re[d" << i << "], Im[d" << i << "], Pwr[d" << i << "], "; + fout << "freq" << fourierData.size()-1 << ", Re[d" << fourierData.size()-1 << "], Im[d" << fourierData.size()-1 << "], Pwr[d" << fourierData.size()-1 << "]" << endl; + + // write data + for (unsigned int j=0; j<data[0].size(); j++) { + for (unsigned int i=0; i<data.size()-1; i++) { + fout << data[i][j] << ", "; + } + fout << data[data.size()-1][j] << endl; + } + fout.close(); + + return 0; +} + +//------------------------------------------------------------------------- +/** + * <p>Groups the histograms before Fourier transform. This is used to group + * detectors. + * + * \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) +{ + // get proper raw run data set + TString runName = *(run.GetRunName()); + PRawRunData *rawRunData = runDataHandler->GetRunData(runName); + if (rawRunData == 0) { + cerr << endl << ">> musrFT_groupHistos **ERROR** Couldn't get raw run data for run '" << runName << "'." << endl; + return 1; + } + + // keep histo list + PIntVector histoList; + for (unsigned int i=0; i<run.GetForwardHistoNoSize(); i++) { + histoList.push_back(run.GetForwardHistoNo(i)); + } + + // check if t0's are found and that #t0 == #histos + PDoubleVector t0; + t0.resize(histoList.size()); + // init t0 vector + for (unsigned int i=0; i<t0.size(); i++) + t0[i] = -1.0; + // 1st: check in the global block + for (unsigned int i=0; i<global->GetT0BinSize(); i++) { + if (i >= t0.size()) { // something is VERY strange + cerr << endl << ">> musrFT_groupHistos **WARNING** found #t0's in GLOBAL block > #histos!"; + cerr << endl << ">> This should NEVER happen. Will ignore these entries."; + cerr << endl << ">> Please check your msr-file!!" << endl; + } else { + t0[i] = global->GetT0Bin(i); + } + } + // 2nd: check in the run block + for (unsigned int i=0; i<run.GetT0BinSize(); i++) { + if (i >= t0.size()) { // something is VERY strange + cerr << endl << ">> musrFT_groupHistos **WARNING** found #t0's in RUN block > #histos!"; + cerr << endl << ">> This should NEVER happen. Will ignore these entries."; + cerr << endl << ">> Please check your msr-file!!" << endl; + } else { + t0[i] = run.GetT0Bin(i); + } + } + // if still some t0's are == -1, estimate t0 + unsigned int idx; + double max; + for (unsigned int i=0; i<t0.size(); i++) { + if (t0[i] == -1.0) { + cout << endl << ">> musrFT_groupHistos **WARNING** try to estimate t0 from maximum in the data set"; + cout << endl << ">> '" << runName << "', histo " << histoList[i] << ". NO warranty this is sensible!"; + idx = 0; + max = rawRunData->GetDataBin(histoList[i])->at(0); + for (unsigned int j=1; j<rawRunData->GetDataBin(histoList[i])->size(); j++) { + if (rawRunData->GetDataBin(histoList[i])->at(j) > max) { + max = rawRunData->GetDataBin(histoList[i])->at(j); + idx = j; + } + } + cout << endl << ">> estimated t0=" << idx << endl; + t0[i] = idx; + } + } + + // group histos + PDoubleVector data = *(rawRunData->GetDataBin(histoList[0])); + for (unsigned int i=1; i<histoList.size(); i++) { + for (unsigned int j=0; j<data.size(); j++) { + if ((j+t0[i]-t0[0] >= 0) && (j+t0[i]-t0[0] < rawRunData->GetDataBin(histoList[i])->size())) { + data[j] += rawRunData->GetDataBin(histoList[i])->at(j); + } + } + } + + rd.rawData.clear(); + rd.rawData = data; + rd.t0 = (int)t0[0]; + + 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; + gettimeofday(&now, 0); + + return ((double)now.tv_sec * 1.0e6 + (double)now.tv_usec)/1.0e3; +} + +//------------------------------------------------------------------------- +/** + * <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 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) { + musrFT_syntax(); + return PMUSR_SUCCESS; + } + + musrFT_startup_param startupParam; + // init startupParam + musrFT_init(startupParam); + + // parse command line options + int status = musrFT_parse_options(argc, argv, startupParam); + if (status != 0) { + int retVal = PMUSR_SUCCESS; + if (status == 2) { + musrFT_syntax(); + retVal = PMUSR_WRONG_STARTUP_SYNTAX; + } + 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; + startup_options.writeExpectedChisq = false; + startup_options.estimateN0 = true; + startup_options.alphaEstimateN0 = 0.0; + TSAXParser *saxParser = new TSAXParser(); + PStartupHandler *startupHandler = new PStartupHandler(); + if (!startupHandler->StartupFileFound()) { + cerr << endl << ">> musrFT **WARNING** couldn't find " << startupHandler->GetStartupFilePath().Data(); + cerr << endl; + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } + } else { + strcpy(startup_path_name, startupHandler->GetStartupFilePath().Data()); + saxParser->ConnectToHandler("PStartupHandler", startupHandler); + //status = saxParser->ParseFile(startup_path_name); + // parsing the file as above seems to lead to problems in certain environments; + // use the parseXmlFile function instead (see PStartupHandler.cpp for the definition) + status = parseXmlFile(saxParser, startup_path_name); + // check for parse errors + if (status) { // error + cerr << endl << ">> musrFT **WARNING** Reading/parsing musrfit_startup.xml failed."; + cerr << endl; + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } + } + } + startupHandler->SetStartupOptions(startup_options); + + // defines the raw time-domain data vector + PPrepFourier data(startupParam.bkg, startupParam.packing); + + // load msr-file(s) + vector<PMsrHandler*> msrHandler; + msrHandler.resize(startupParam.msrFln.size()); + for (unsigned int i=0; i<startupParam.msrFln.size(); i++) { + msrHandler[i] = new PMsrHandler(startupParam.msrFln[i].Data(), startupHandler->GetStartupOptions(), true); + status = msrHandler[i]->ReadMsrFile(); + if (status != PMUSR_SUCCESS) { + switch (status) { + case PMUSR_MSR_FILE_NOT_FOUND: + cout << endl << ">> musrFT **ERROR** couldn't find " << startupParam.msrFln[i] << endl << endl; + break; + case PMUSR_MSR_SYNTAX_ERROR: + cout << endl << ">> musrFT **SYNTAX ERROR** in file " << startupParam.msrFln[i] << ", full stop here." << endl << endl; + break; + default: + cout << endl << ">> musrFT **UNKOWN ERROR** when trying to read the msr-file" << endl << endl; + break; + } + return status; + } + } + + vector<PRunDataHandler*> runDataHandler; + runDataHandler.resize(startupParam.msrFln.size()+startupParam.dataFln.size()); // resize to the total number of run data provided + // load data-file(s) related to msr-file + for (unsigned int i=0; i<msrHandler.size(); i++) { + // create run data handler + if (startupHandler) + runDataHandler[i] = new PRunDataHandler(msrHandler[i], startupHandler->GetDataPathList()); + else + runDataHandler[i] = new PRunDataHandler(msrHandler[i]); + } + + // load data-file(s) provided directly + for (unsigned int i=msrHandler.size(); i<msrHandler.size()+startupParam.dataFln.size(); i++) { + // create run data handler + if (startupHandler) + runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList()); + else + runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()]); + } + + // read all the data files + for (unsigned int i=0; i<runDataHandler.size(); i++) { + runDataHandler[i]->ReadData(); + + if (!runDataHandler[i]->IsAllDataAvailable()) { + if (i < msrHandler.size()) { + cerr << endl << ">> musrFT **ERROR** couldn't read data from msr-file '" << startupParam.msrFln[i] << "'." << endl; + } else { + cerr << endl << ">> musrFT **ERROR** couldn't read data-file '" << startupParam.dataFln[i] << "'." << endl; + } + return PMUSR_DATA_FILE_READ_ERROR; + } + + // dig out all the necessary time domain data + PRawRunData *rawRunData = runDataHandler[i]->GetRunData(); + if (rawRunData == 0) { + if (i < msrHandler.size()) { + cerr << endl << ">> musrFT **ERROR** couldn't obtain the raw run data set from msr-file " << startupParam.msrFln[i] << endl; + } else { + cerr << endl << ">> musrFT **ERROR** couldn't obtain the raw run data set for " << startupParam.dataFln[i-msrHandler.size()] << endl; + } + return PMUSR_DATA_FILE_READ_ERROR; + } + + // first check of histo list makes sense + if (i >= msrHandler.size()) { // only check if originating from data-files (not msr-files) + for (unsigned int j=0; j<startupParam.histo.size(); j++) { + if ((unsigned int)startupParam.histo[j] > rawRunData->GetNoOfHistos()) { + cerr << endl << ">> musrFT **ERROR** found histo no " << startupParam.histo[j] << " > # of histo in the file ("; + cerr << startupParam.dataFln[i] << " // # histo: " << rawRunData->GetNoOfHistos() << ")." << endl; + return PMUSR_DATA_FILE_READ_ERROR; + } + } + if (startupParam.histo.size() == 0) { // no histo list given + // set histo list to ALL available histos for the data file + for (unsigned int j=0; j<rawRunData->GetNoOfHistos(); j++) + startupParam.histo.push_back(j+1); + } + } + + musrFT_data rd; + TString str(""), fln(""); + 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) { + cerr << endl << ">> musrFT **ERROR** couldn't obtain PLOT block from msr-handler." << endl; + return PMUSR_DATA_FILE_READ_ERROR; + } + // keep RUN block(s) info + PMsrRunList *runs = msrHandler[i]->GetMsrRunList(); + if (runs == 0) { + cerr << endl << ">> musrFT **ERROR** couldn't obtain RUN block(s) from msr-handler." << endl; + return PMUSR_DATA_FILE_READ_ERROR; + } + // keep GLOBAL block info + PMsrGlobalBlock *global = msrHandler[i]->GetMsrGlobal(); + if (global == 0) { + cerr << endl << ">> musrFT **ERROR** couldn't obtain GLOBAL block from msr-handler." << endl; + return PMUSR_DATA_FILE_READ_ERROR; + } + // keep FOURIER block info + PMsrFourierStructure *fourierBlock = msrHandler[i]->GetMsrFourierList(); + if (fourierBlock == 0) { + cerr << endl << ">> msrFT **WARNING** couldn't obtain FOURIER block from msr-handler." << endl; + return PMUSR_DATA_FILE_READ_ERROR; + } else { // filter out all necessary info + if (fourierBlock->fFourierBlockPresent) { + // get units + unitTag = fourierBlock->fUnits; + // get fourier power + if (startupParam.fourierPower == -1) { // no Fourier power given from the command line, hence check FOURIER block + if (fourierBlock->fFourierPower > 1) + startupParam.fourierPower = fourierBlock->fFourierPower; + } + // get apodization tag + apodTag = fourierBlock->fApodization; + // 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]; + } + // 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; + } + } + } + + // get the run information from the msr-file PLOT block 'runs' + PIntVector runList = plot->at(0).fRuns; + + // loop over all runs listed in the msr-file PLOT block + for (unsigned int j=0; j<runList.size(); j++) { + + // keep forward histo list + PIntVector histoList; + for (unsigned int k=0; k<runs->at(runList[j]-1).GetForwardHistoNoSize(); k++) { + histoList.push_back(runs->at(runList[j]-1).GetForwardHistoNo(k)); + } + + // handle meta information + fln = runDataHandler[i]->GetRunPathName(); + musrFT_getMetaInfo(fln, rawRunData, str); + TString hh(""); + hh = TString::Format("h%d", histoList[0]); + for (unsigned int k=1; k<histoList.size(); k++) + hh += TString::Format("/%d", histoList[k]); + hh += ":"; + rd.info = hh; + rd.info += str; + + // handle time resolution + rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us) + + // handle time range + // take it from msr-file PLOT block 'range' if not overwritten from the command line + if ((startupParam.timeRange[0] != -1) && (startupParam.timeRange[1] != -1)) { + rd.timeRange[0] = startupParam.timeRange[0]; + rd.timeRange[1] = startupParam.timeRange[1]; + } else { + 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) + // group forward histos + if (musrFT_groupHistos(runDataHandler[i], global, runs->at(runList[j]-1), rd)) { + return PMUSR_DATA_FILE_READ_ERROR; + } + // keep data set + data.AddData(rd); + + // get packing + Int_t pack = 1; + if (global->GetPacking() != -1) { + pack = global->GetPacking(); + } + if (runs->at(runList[j]-1).GetPacking() != -1) { + pack = runs->at(runList[j]-1).GetPacking(); + } + if (startupParam.packing > 1) + pack = startupParam.packing; + data.SetPacking(pack); + } + } else { // obtain info from command line options for direct data-file read + musrFT_getMetaInfo(startupParam.dataFln[i-msrHandler.size()], rawRunData, str); + for (unsigned int j=0; j<startupParam.histo.size(); j++) { + idx = startupParam.histo[j]; + + // handle meta information + rd.info = TString::Format("h%d:", idx); + rd.info += str; + + // handle time resolution + rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us) + + // handle time range + rd.timeRange[0] = startupParam.timeRange[0]; // in (us) + rd.timeRange[1] = startupParam.timeRange[1]; // in (us) + + // handle data set + rd.rawData.clear(); + rd.rawData = *(rawRunData->GetDataBin(idx)); + + // handle t0's + rd.t0 = -1; + if (startupParam.t0.size() == 1) + rd.t0 = startupParam.t0[0]; + else if (j < startupParam.t0.size()) + rd.t0 = startupParam.t0[j]; + if (rd.t0 == -1) { // no t0 given, try to estimate it + musrFT_estimateT0(rd); + } + + data.AddData(rd); + } + } + } + + // 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(); + + // do the time domain filtering now + data.DoFiltering(); + + // do lifetime correction + if (startupParam.lifetimecorrection != 0.0) + data.DoLifeTimeCorrection(startupParam.lifetimecorrection); + + // do packing + data.DoPacking(); + + // get all the corrected data + vector<TH1F*> histo = data.GetData(); + + // prepare Fourier + if (startupParam.fourierUnits.BeginsWith("gauss", TString::kIgnoreCase)) + unitTag = FOURIER_UNIT_GAUSS; + else if (startupParam.fourierUnits.BeginsWith("tesla", TString::kIgnoreCase)) + unitTag = FOURIER_UNIT_TESLA; + else if (startupParam.fourierUnits.BeginsWith("mhz", TString::kIgnoreCase)) + unitTag = FOURIER_UNIT_FREQ; + else if (startupParam.fourierUnits.BeginsWith("mc/s", TString::kIgnoreCase)) + unitTag = FOURIER_UNIT_CYCLES; + else if (startupParam.fourierUnits.BeginsWith("??", TString::kIgnoreCase) && (unitTag == FOURIER_UNIT_NOT_GIVEN)) + unitTag = FOURIER_UNIT_FREQ; + + vector<PFourier*> fourier; + fourier.resize(histo.size()); + for (unsigned int i=0; i<fourier.size(); i++) { + fourier[i] = new PFourier(histo[i], unitTag, 0.0, 0.0, true, startupParam.fourierPower); + } + + // Fourier transform data + if (startupParam.apodization.BeginsWith("weak", TString::kIgnoreCase)) + apodTag = F_APODIZATION_WEAK; + else if (startupParam.apodization.BeginsWith("medium", TString::kIgnoreCase)) + apodTag = F_APODIZATION_MEDIUM; + else if (startupParam.apodization.BeginsWith("strong", TString::kIgnoreCase)) + apodTag = F_APODIZATION_STRONG; + + double start = millitime(); + for (unsigned int i=0; i<fourier.size(); i++) { + fourier[i]->Transform(apodTag); + } + 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 { // do Canvas + + // 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 + 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; + + for (unsigned int i=0; i<msrHandler.size(); i++) + if (msrHandler[i]) + delete msrHandler[i]; + msrHandler.clear(); + + for (unsigned int i=0; i<runDataHandler.size(); i++) + if (runDataHandler[i]) + delete runDataHandler[i]; + runDataHandler.clear(); + + if (histo.size() > 0) { + for (unsigned int i=0; i<histo.size(); i++) + delete histo[i]; + histo.clear(); + } + if (fourier.size() > 0) { + for (unsigned int i=0; i<fourier.size(); i++) + delete fourier[i]; + fourier.clear(); + } + + return PMUSR_SUCCESS; +} diff --git a/src/musrview.cpp b/src/musrview.cpp index e9cb2558..17d3d978 100644 --- a/src/musrview.cpp +++ b/src/musrview.cpp @@ -67,6 +67,8 @@ void musrview_syntax() cout << endl << " eps, pdf, gif, jpg, png, svg, xpm, root"; cout << endl << " example: musrview 3310.msr --png, will produce a files 3310_X.png"; cout << endl << " where 'X' stands for the plot number (starting form 0)"; + cout << endl << " --ascii: "; + cout << endl << " will produce an ascii dump of the data and fit as plotted."; cout << endl << " --timeout <timeout>: <timeout> given in seconds after which musrview terminates."; cout << endl << " If <timeout> <= 0, no timeout will take place. Default <timeout> is 0."; cout << endl; @@ -100,6 +102,7 @@ int main(int argc, char *argv[]) bool success = true; char fileName[128]; bool graphicsOutput = false; + bool asciiOutput = false; char graphicsExtension[128]; int timeout = 0; @@ -135,6 +138,8 @@ int main(int argc, char *argv[]) graphicsOutput = true; strcpy(graphicsExtension, argv[i]+2); + } else if (!strcmp(argv[i], "--ascii")) { + asciiOutput = true; } else if (!strcmp(argv[i], "--timeout")) { if (i+1 < argc) { TString str(argv[i+1]); @@ -280,7 +285,7 @@ int main(int argc, char *argv[]) if (success) { // generate Root application needed for PMusrCanvas - if (graphicsOutput) { + if (graphicsOutput || asciiOutput) { argv[argc] = (char*)malloc(16*sizeof(char)); strcpy(argv[argc], "-b"); argc++; @@ -299,10 +304,10 @@ int main(int argc, char *argv[]) startupHandler->GetFourierDefaults(), startupHandler->GetMarkerList(), startupHandler->GetColorList(), - graphicsOutput); + graphicsOutput||asciiOutput); else musrCanvas = new PMusrCanvas(i, msrHandler->GetMsrTitle()->Data(), - 10+i*100, 10+i*100, 800, 600, graphicsOutput); + 10+i*100, 10+i*100, 800, 600, graphicsOutput||asciiOutput); if (!musrCanvas->IsValid()) { cerr << endl << ">> musrview **SEVERE ERROR** Couldn't invoke all necessary objects, will quit."; @@ -334,6 +339,12 @@ int main(int argc, char *argv[]) musrCanvas->SaveGraphicsAndQuit(fileName, graphicsExtension); } + if (asciiOutput) { + // save data in batch mode + musrCanvas->SaveDataAscii(); + musrCanvas->Done(0); + } + // keep musrCanvas objects canvasVector.push_back(musrCanvas); } @@ -349,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();