From c6900030c7ff8dcb0e41d020657197bf73914139 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Mon, 2 Feb 2015 15:52:51 +0100 Subject: [PATCH 1/5] first partial implementation for a standalone Fourier transform calculator/plotter. Mainly meant for HAL-9500. Much more to come ... --- ChangeLog | 4 + src/Makefile.am | 3 +- src/classes/Makefile.am | 2 + src/classes/PPrepFourier.cpp | 389 +++++++++++++ src/classes/PRunDataHandler.cpp | 21 + src/include/PFourier.h | 4 +- src/include/PMusr.h | 3 +- src/include/PPrepFourier.h | 86 +++ src/include/PRunDataHandler.h | 3 +- src/musrFT.cpp | 967 ++++++++++++++++++++++++++++++++ 10 files changed, 1477 insertions(+), 5 deletions(-) create mode 100644 src/classes/PPrepFourier.cpp create mode 100644 src/include/PPrepFourier.h create mode 100644 src/musrFT.cpp diff --git a/ChangeLog b/ChangeLog index bcc8ab9b..4f92c0c7 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,6 +4,10 @@ changes since 0.13.0 =================================== +NEW 2015-02-02 first partial implementation of a standalone Fourier transform/plotter: + musrFT. A LOT of the functionality is still missing! Initially it is + meant to be used for HAL-9500, i.e. Fourier transform WITHOUT lifetime + correction. NEW 2014-12-18 first implementation of a GLOBAL block which allows to shorten a typical msr-file. Duplicate entries from the RUN blocks can be added here. Furthermore, the 'lifetimecorrection' flag is diff --git a/src/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 51283653..0be0160e 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -12,6 +12,7 @@ h_sources = \ ../include/PMusrCanvas.h \ ../include/PMusr.h \ ../include/PMusrT0.h \ + ../include/PPrepFourier.h \ ../include/PRunAsymmetry.h \ ../include/PRunBase.h \ ../include/PRunDataHandler.h \ @@ -52,6 +53,7 @@ cpp_sources = \ PMusrCanvas.cpp \ PMusr.cpp \ PMusrT0.cpp \ + PPrepFourier.cpp \ PRunAsymmetry.cpp \ PRunBase.cpp \ PRunDataHandler.cpp \ diff --git a/src/classes/PPrepFourier.cpp b/src/classes/PPrepFourier.cpp new file mode 100644 index 00000000..a333d625 --- /dev/null +++ b/src/classes/PPrepFourier.cpp @@ -0,0 +1,389 @@ +/*************************************************************************** + + 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 +//-------------------------------------------------------------------------- +/** + *

+ * + * \param bkgRange + */ +void PPrepFourier::SetBkgRange(const Int_t *bkgRange) +{ + cout << endl << "debug> bkgRange: " << bkgRange[0] << ", " << bkgRange[1] << endl; + + int err=0; + if (bkgRange[0] >= -1) { + fBkgRange[0] = bkgRange[0]; + } 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 +//-------------------------------------------------------------------------- +/** + *

+ * + * \param packing + */ +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 +//-------------------------------------------------------------------------- +/** + *

+ * + * \param data + */ +void PPrepFourier::AddData(musrFT_data &data) +{ + fRawData.push_back(data); +} + +//-------------------------------------------------------------------------- +// DoBkgCorrection +//-------------------------------------------------------------------------- +/** + *

+ * + */ +void PPrepFourier::DoBkgCorrection() +{ + cout << endl << "debug> DoBkgCorrection ..."; + + // if no bkg-range is given, nothing needs to be done + if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) { + cout << endl << "debug> no background range given ..."; + return; + } + + // 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; + } + } + + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } + + Double_t bkg=0.0; + for (unsigned int i=0; i histo: " << i << ", bkg=" << bkg; + + // correct data + for (unsigned int j=0; j + * + */ +void PPrepFourier::DoPacking() +{ + cout << endl << "debug> DoPacking ..."; + + if (fPacking == 1) // nothing to be done + return; + + PDoubleVector tmpData; + Double_t dval = 0.0; + for (unsigned int i=0; i histo " << i+1 << ": packed data: "; + for (unsigned int j=0; j<15; j++) + cout << fData[i][j] << ", "; + } +} + +//-------------------------------------------------------------------------- +// DoFiltering +//-------------------------------------------------------------------------- +/** + *

+ * + */ +void PPrepFourier::DoFiltering() +{ + cout << endl << "debug> DoFiltering not yet implemented ..."; +} + +//-------------------------------------------------------------------------- +// GetInfo +//-------------------------------------------------------------------------- +/** + *

+ * + * \param idx + */ +TString PPrepFourier::GetInfo(const UInt_t idx) +{ + TString info(""); + + if (idx < fRawData.size()) + info = fRawData[idx].info; + + return info; +} + +//-------------------------------------------------------------------------- +// GetData +//-------------------------------------------------------------------------- +/** + *

+ * + */ +vector PPrepFourier::GetData() +{ + vector data; + data.resize(fData.size()); + + 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 + */ +TH1F *PPrepFourier::GetData(const UInt_t idx) +{ + if (idx > fData.size()) + 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()); + for (unsigned int i=0; iChecks 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/include/PFourier.h b/src/include/PFourier.h index e8355290..acb7a5b2 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 * @@ -68,7 +68,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/PMusr.h b/src/include/PMusr.h index f9826fbd..b4fe1034 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 diff --git a/src/include/PPrepFourier.h b/src/include/PPrepFourier.h new file mode 100644 index 00000000..852b874e --- /dev/null +++ b/src/include/PPrepFourier.h @@ -0,0 +1,86 @@ +/*************************************************************************** + + 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" + +//---------------------------------------------------------------------------- +/** + *

+ */ +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(); + + 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..9a13c261 --- /dev/null +++ b/src/musrFT.cpp @@ -0,0 +1,967 @@ +/*************************************************************************** + + 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 "git-revision.h" +#include "PMusr.h" +#include "PStartupHandler.h" +#include "PMsrHandler.h" +#include "PRunDataHandler.h" +#include "PPrepFourier.h" +#include "PFourier.h" + +//---------------------------------------------------------------------------- +/** + *

+ */ +typedef struct { + vector msrFln; + vector dataFln; + vector dataFileFormat; + TString graphicFormat; + TString dumpFln; + int bkg[2]; + TString fourierOpt; + TString apotization; + int fourierPower; + TString fourierUnits; + double initialPhase; + double fourierRange[2]; + double timeRange[2]; + vector histo; + bool showAverage; + vector t0; + int packing; + TString title; +} musrFT_startup_param; + +//------------------------------------------------------------------------- +/** + *

+ */ +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', 'power', 'imag', 'real+imag', of 'phase'."; + cout << endl << " If this is not defined (neither on the command line nor in the musrFT_startup.xml),"; + cout << endl << " default will be 'real'."; + cout << endl << " -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 inital 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 << " -t, --title : give a global title for the plot."; + cout << endl << endl; +} + +//------------------------------------------------------------------------- +/** + * <p>initialize startup parameters. + * + * \param startupParam + */ +void musrFT_init(musrFT_startup_param &startupParam) +{ + startupParam.graphicFormat = TString(""); + startupParam.dumpFln = TString(""); + startupParam.bkg[0] = -1; + startupParam.bkg[1] = -1; + startupParam.fourierOpt = TString("real"); + startupParam.apotization = TString("none"); + startupParam.fourierPower = -1; + startupParam.fourierUnits = TString("MHz"); + 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(""); +} + +//------------------------------------------------------------------------- +/** + * <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 + * \param argv + * \param startupParam + */ +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.apotization = 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("-t ") || 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("-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; +} + +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(); +} + +//------------------------------------------------------------------------- +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 << "debug> estimated t0=" << idx << endl; + rd.t0 = idx; +} + +//------------------------------------------------------------------------- +void musrFT_cleanup(TH1F *h) +/** + * <p> + * + * \param h histogram to be deleted + */ +{ + if (h) { + delete h; + h = 0; + } +} + +//------------------------------------------------------------------------- +/** + * <p> + * + * \param data + */ +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; +} + +//--------------------------------------------------- +double millitime() +{ + struct timeval now; + gettimeofday(&now, 0); + + return ((double)now.tv_sec * 1.0e6 + (double)now.tv_usec)/1.0e3; +} + +//------------------------------------------------------------------------- +/** + * <p> + * + * \param argc + * \param argv + */ +int main(int argc, char *argv[]) +{ + // only program name alone + if (argc == 1) { + musrFT_syntax(); + return PMUSR_SUCCESS; + } + + musrFT_startup_param startupParam; + // init startupParam + musrFT_init(startupParam); + + // parse 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; + } + + // 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) + cout << endl << "debug> before reading msr-files ..."; + 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; + } + } + + cout << endl << "debug> before reading data-files ..."; + 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++) { + cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()]; + // create run data handler + if (startupHandler) + runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList()); + else + runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()]); + + // read the data files + runDataHandler[i]->ReadData(); + + if (!runDataHandler[i]->IsAllDataAvailable()) { + cerr << endl << ">> musrFT **ERROR** couldn't read data-file '" << startupParam.dataFln[i-msrHandler.size()] << "'."; + return PMUSR_DATA_FILE_READ_ERROR; + } + + // dig out all the necessary time domain data + PRawRunData *rawRunData = runDataHandler[i]->GetRunData(); + if (rawRunData == 0) { + 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 + 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-msrHandler.size()] << " // # histo: " << rawRunData->GetNoOfHistos() << ")." << endl; + return PMUSR_DATA_FILE_READ_ERROR; + } + } + + musrFT_data rd; + TString str(""); + unsigned int idx=0; + musrFT_getMetaInfo(startupParam.dataFln[i-msrHandler.size()], rawRunData, str); + for (unsigned int j=0; j<startupParam.histo.size(); j++) { + idx = startupParam.histo[j]; + // meta information + rd.info = TString::Format("h%d:", idx); + rd.info += str; +/* + cout << endl << "debug> ++++++++"; + cout << endl << "debug> info: " << rd.info; +*/ + // time resolution + rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us) +// cout << endl << "debug> time resolution: " << rd.timeResolution << " (us)"; + // time range + rd.timeRange[0] = startupParam.timeRange[0]; // in (us) + rd.timeRange[1] = startupParam.timeRange[1]; // in (us) +// cout << endl << "debug> time range: " << rd.timeRange[0] << ".." << rd.timeRange[1] << " (us)"; + // data set + rd.rawData.clear(); + rd.rawData = *(rawRunData->GetDataBin(idx)); + // check for potentially given 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]; +// cout << endl << "debug> t0 = " << rd.t0; + if (rd.t0 == -1) { // no t0 given, try to estimate it + musrFT_estimateT0(rd); + } +/* + cout << endl << "debug> rd.rawData.size()=" << rd.rawData.size(); + cout << endl << "debug> rd.rawData: h" << idx << ": "; +*/ + for (unsigned int k=rd.t0; k<rd.t0+15; k++) + cout << rd.rawData[k] << ", "; + data.AddData(rd); + } + } + +/* + cout << endl << "debug> +++++++++"; + cout << endl << "debug> data.GetNoOfData()=" << data.GetNoOfData(); +*/ + + // calculate background levels and subtract them from the data + data.DoBkgCorrection(); + + // do the time domain filtering now + data.DoFiltering(); + + // do packing + data.DoPacking(); + + // get all the corrected data + vector<TH1F*> histo = data.GetData(); +/* + if (histo.size() == 0) + cout << endl << "debug> histo.size()==0!! :-(" << endl; + for (unsigned int i=0; i<histo.size(); i++) { + cout << endl << "debug> histo: time resolution " << histo[i]->GetBinWidth(1) << " (ns)."; + cout << endl << "debug> histo[" << i << "]->GetNbinsX()=" << histo[i]->GetNbinsX(); + cout << endl << "debug> histo data: "; + for (unsigned int j=1; j<=15; j++) + cout << histo[i]->GetBinContent(j) << ", "; + cout << endl; + } +*/ + // prepare Fourier + Int_t unitTag = FOURIER_UNIT_NOT_GIVEN; + 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; + + 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 + Int_t apodTag = F_APODIZATION_NONE; + if (startupParam.apotization.BeginsWith("weak", TString::kIgnoreCase)) + apodTag = F_APODIZATION_WEAK; + else if (startupParam.apotization.BeginsWith("medium", TString::kIgnoreCase)) + apodTag = F_APODIZATION_MEDIUM; + else if (startupParam.apotization.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; + + // 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 { + + // if Fourier graphical export is whished, switch to batch mode + + // plot the Fourier transform + + } + + + // cleanup + 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; +} From d4ba2f6d81150631a430eb50aa914b188cfcdce4 Mon Sep 17 00:00:00 2001 From: Suter Andreas <andreas.suter@psi.ch> Date: Wed, 11 Feb 2015 15:35:14 +0100 Subject: [PATCH 2/5] more work towards a standalone Fourier (mainly for HAL-9500). Add msr-file support. Still some severe flaws. --- src/classes/PMsrHandler.cpp | 44 +++-- src/classes/PPrepFourier.cpp | 80 +++++++- src/include/PMusr.h | 2 +- src/include/PPrepFourier.h | 1 + src/musrFT.cpp | 364 +++++++++++++++++++++++++++-------- 5 files changed, 396 insertions(+), 95 deletions(-) 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; i<param.GetMap()->size(); 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; i<param.GetMap()->size(); 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/PPrepFourier.cpp b/src/classes/PPrepFourier.cpp index a333d625..30750f8f 100644 --- a/src/classes/PPrepFourier.cpp +++ b/src/classes/PPrepFourier.cpp @@ -155,6 +155,11 @@ void PPrepFourier::DoBkgCorrection() { cout << endl << "debug> DoBkgCorrection ..."; + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } + // if no bkg-range is given, nothing needs to be done if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) { cout << endl << "debug> no background range given ..."; @@ -169,11 +174,6 @@ void PPrepFourier::DoBkgCorrection() } } - // make sure fData are already present, and if not create the necessary data sets - if (fData.size() != fRawData.size()) { - InitData(); - } - Double_t bkg=0.0; for (unsigned int i=0; i<fRawData.size(); i++) { // calculate the bkg for the given range @@ -198,7 +198,12 @@ void PPrepFourier::DoBkgCorrection() */ void PPrepFourier::DoPacking() { - cout << endl << "debug> DoPacking ..."; + cout << endl << "debug> DoPacking : pack=" << fPacking << " ..."; + + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } if (fPacking == 1) // nothing to be done return; @@ -236,6 +241,53 @@ void PPrepFourier::DoPacking() void PPrepFourier::DoFiltering() { cout << endl << "debug> DoFiltering not yet implemented ..."; + + // make sure fData are already present, and if not create the necessary data sets + if (fData.size() != fRawData.size()) { + InitData(); + } + +} + +//-------------------------------------------------------------------------- +// DoLifeTimeCorrection +//-------------------------------------------------------------------------- +/** + * <p> + * + * \param fudge + */ +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; i<fData.size(); i++) { + scale = fRawData[i].timeResolution / PMUON_LIFETIME; + for (unsigned int j=0; j<fData[i].size(); j++) { + fData[i][j] *= exp(j*scale); + } + } + + // calc N0 + Double_t dval; + Double_t N0; + for (unsigned int i=0; i<fData.size(); i++) { + dval = 0.0; + for (unsigned int j=0; j<fData[i].size(); j++) { + dval += fData[i][j]; + } + N0 = dval/fData[i].size(); + N0 *= fudge; + for (unsigned int j=0; j<fData[i].size(); j++) { + fData[i][j] -= N0; + fData[i][j] /= N0; + } + } } //-------------------------------------------------------------------------- @@ -268,6 +320,10 @@ vector<TH1F*> PPrepFourier::GetData() vector<TH1F*> 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; @@ -330,7 +386,10 @@ vector<TH1F*> PPrepFourier::GetData() */ TH1F *PPrepFourier::GetData(const UInt_t idx) { - if (idx > fData.size()) + 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); @@ -381,8 +440,13 @@ TH1F *PPrepFourier::GetData(const UInt_t idx) void PPrepFourier::InitData() { fData.resize(fRawData.size()); + unsigned int t0; for (unsigned int i=0; i<fRawData.size(); i++) { - for (unsigned int j=fRawData[i].t0; j<fRawData[i].rawData.size(); j++) { + if (fRawData[i].t0 >= 0) + t0 = fRawData[i].t0; + else + t0 = 0; + for (unsigned int j=t0; j<fRawData[i].rawData.size(); j++) { fData[i].push_back(fRawData[i].rawData[j]); } } diff --git a/src/include/PMusr.h b/src/include/PMusr.h index b4fe1034..f083b13b 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -704,7 +704,7 @@ typedef vector<PMsrRunBlock> 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/PPrepFourier.h b/src/include/PPrepFourier.h index 852b874e..735a23b7 100644 --- a/src/include/PPrepFourier.h +++ b/src/include/PPrepFourier.h @@ -67,6 +67,7 @@ class PPrepFourier { void DoBkgCorrection(); void DoPacking(); void DoFiltering(); + void DoLifeTimeCorrection(Double_t fudge); TString GetInfo(const UInt_t idx); UInt_t GetNoOfData() { return fRawData.size(); } diff --git a/src/musrFT.cpp b/src/musrFT.cpp index 9a13c261..3551eedb 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -74,6 +74,7 @@ typedef struct { vector<int> t0; int packing; TString title; + double lifetimecorrection; //< is == 0.0 for NO life time correction, otherwise it holds the fudge factor } musrFT_startup_param; //------------------------------------------------------------------------- @@ -100,7 +101,7 @@ void musrFT_syntax() cout << endl << " Supported graphic-format-extension: eps, pdf, gif, jpg, png, svg, xpm, root"; cout << endl << " --dump <fln> : 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 <fln>."; - cout << endl << " --filter : filter and filter-specific-information -- TO BE WRITTEN YET."; + cout << endl << " --filter : filter and filter-specific-information -- ***TO BE WRITTEN YET***."; cout << endl << " -b, --background <start> <end>: background interval used to estimate the backround to be"; cout << endl << " subtracted before the Fourier transform. <start>, <end> to be given in bins."; cout << endl << " -fo, --fourier-option <fopt>: <fopt> can be 'real', 'power', 'imag', 'real+imag', of 'phase'."; @@ -128,12 +129,17 @@ void musrFT_syntax() 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 <list> : 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 << " --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 <N> : if <N> (an integer), the time domain data will first be packed/rebinned by <N>."; cout << endl << " -t, --title <title> : give a global title for the plot."; + cout << endl << " --create-msr-file <fln> : creates a msr-file based on the command line options"; + cout << endl << " provided. This will help on the way to a full fitting model."; + cout << endl << " ***TO BE WRITTEN YET.***"; + cout << endl << " -lc, --lifetimecorrection <fudge>: try to eliminate muon life time decay. Only makes sense for low"; + cout << endl << " transverse fields. <fudge> is a tweaking factor and should be kept around 1.0."; cout << endl << endl; } @@ -152,7 +158,7 @@ void musrFT_init(musrFT_startup_param &startupParam) startupParam.fourierOpt = TString("real"); startupParam.apotization = TString("none"); startupParam.fourierPower = -1; - startupParam.fourierUnits = TString("MHz"); + startupParam.fourierUnits = TString("??"); startupParam.initialPhase = 0.0; startupParam.fourierRange[0] = -1.0; startupParam.fourierRange[1] = -1.0; @@ -161,6 +167,7 @@ void musrFT_init(musrFT_startup_param &startupParam) startupParam.showAverage = false; startupParam.packing = 1; startupParam.title = TString(""); + startupParam.lifetimecorrection = 0.0; } //------------------------------------------------------------------------- @@ -461,6 +468,18 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa return 2; } startupParam.packing = pack.Atoi(); + } 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 an double '" << fudge << "'." << endl; + return 2; + } + startupParam.lifetimecorrection = fudge.Atof(); } else if (tstr.BeginsWith("-df") || tstr.BeginsWith("--data-file")) { while (++i < argc) { if (argv[i][0] == '-') { @@ -520,6 +539,10 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa return 0; } +//---------------------------------------------------------------------------------------- +/** + * + */ void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &metaInfo) { double dval; @@ -572,7 +595,7 @@ void musrFT_estimateT0(musrFT_data &rd) idx = (int)i; } } - cout << endl << "debug> estimated t0=" << idx << endl; + cout << endl << ">> musrFT_estimateT0: estimated t0=" << idx << endl; rd.t0 = idx; } @@ -668,6 +691,94 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d return 0; } +//--------------------------------------------------- +/** + * <p> + * + * \param runDataHandler + * \param global + * \param run + * \param rd + */ +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; +} + //--------------------------------------------------- double millitime() { @@ -686,6 +797,9 @@ double millitime() */ int main(int argc, char *argv[]) { + Int_t unitTag = FOURIER_UNIT_NOT_GIVEN; + Int_t apodTag = F_APODIZATION_NONE; + // only program name alone if (argc == 1) { musrFT_syntax(); @@ -696,7 +810,7 @@ int main(int argc, char *argv[]) // init startupParam musrFT_init(startupParam); - // parse options + // parse command line options int status = musrFT_parse_options(argc, argv, startupParam); if (status != 0) { int retVal = PMUSR_SUCCESS; @@ -755,7 +869,6 @@ int main(int argc, char *argv[]) PPrepFourier data(startupParam.bkg, startupParam.packing); // load msr-file(s) - cout << endl << "debug> before reading msr-files ..."; vector<PMsrHandler*> msrHandler; msrHandler.resize(startupParam.msrFln.size()); for (unsigned int i=0; i<startupParam.msrFln.size(); i++) { @@ -777,7 +890,6 @@ int main(int argc, char *argv[]) } } - cout << endl << "debug> before reading data-files ..."; 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 @@ -791,110 +903,209 @@ int main(int argc, char *argv[]) // load data-file(s) provided directly for (unsigned int i=msrHandler.size(); i<msrHandler.size()+startupParam.dataFln.size(); i++) { - cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()]; +// cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()]; // create run data handler if (startupHandler) runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList()); else runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()]); + } - // read the data files + // read all the data files + for (unsigned int i=0; i<runDataHandler.size(); i++) { runDataHandler[i]->ReadData(); if (!runDataHandler[i]->IsAllDataAvailable()) { - cerr << endl << ">> musrFT **ERROR** couldn't read data-file '" << startupParam.dataFln[i-msrHandler.size()] << "'."; + 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) { - cerr << endl << ">> musrFT **ERROR** couldn't obtain the raw run data set for " << startupParam.dataFln[i-msrHandler.size()] << endl; + 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 - 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-msrHandler.size()] << " // # histo: " << rawRunData->GetNoOfHistos() << ")." << endl; - return PMUSR_DATA_FILE_READ_ERROR; + 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(""); + TString str(""), fln(""); unsigned int idx=0; - musrFT_getMetaInfo(startupParam.dataFln[i-msrHandler.size()], rawRunData, str); - for (unsigned int j=0; j<startupParam.histo.size(); j++) { - idx = startupParam.histo[j]; - // meta information - rd.info = TString::Format("h%d:", idx); - rd.info += str; -/* - cout << endl << "debug> ++++++++"; - cout << endl << "debug> info: " << rd.info; -*/ - // time resolution - rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us) -// cout << endl << "debug> time resolution: " << rd.timeResolution << " (us)"; - // time range - rd.timeRange[0] = startupParam.timeRange[0]; // in (us) - rd.timeRange[1] = startupParam.timeRange[1]; // in (us) -// cout << endl << "debug> time range: " << rd.timeRange[0] << ".." << rd.timeRange[1] << " (us)"; - // data set - rd.rawData.clear(); - rd.rawData = *(rawRunData->GetDataBin(idx)); - // check for potentially given 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]; -// cout << endl << "debug> t0 = " << rd.t0; - if (rd.t0 == -1) { // no t0 given, try to estimate it - musrFT_estimateT0(rd); + // get meta info, time resolution, time range, raw data sets + if (i < msrHandler.size()) { // obtain info from msr-files + // 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 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 { + 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); } -/* - cout << endl << "debug> rd.rawData.size()=" << rd.rawData.size(); - cout << endl << "debug> rd.rawData: h" << idx << ": "; -*/ - for (unsigned int k=rd.t0; k<rd.t0+15; k++) - cout << rd.rawData[k] << ", "; - data.AddData(rd); } } -/* - cout << endl << "debug> +++++++++"; - cout << endl << "debug> data.GetNoOfData()=" << data.GetNoOfData(); -*/ - // 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(); -/* - if (histo.size() == 0) - cout << endl << "debug> histo.size()==0!! :-(" << endl; - for (unsigned int i=0; i<histo.size(); i++) { - cout << endl << "debug> histo: time resolution " << histo[i]->GetBinWidth(1) << " (ns)."; - cout << endl << "debug> histo[" << i << "]->GetNbinsX()=" << histo[i]->GetNbinsX(); - cout << endl << "debug> histo data: "; - for (unsigned int j=1; j<=15; j++) - cout << histo[i]->GetBinContent(j) << ", "; - cout << endl; - } -*/ + // prepare Fourier - Int_t unitTag = FOURIER_UNIT_NOT_GIVEN; if (startupParam.fourierUnits.BeginsWith("gauss", TString::kIgnoreCase)) unitTag = FOURIER_UNIT_GAUSS; else if (startupParam.fourierUnits.BeginsWith("tesla", TString::kIgnoreCase)) @@ -903,6 +1114,8 @@ int main(int argc, char *argv[]) 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()); @@ -911,7 +1124,6 @@ int main(int argc, char *argv[]) } // Fourier transform data - Int_t apodTag = F_APODIZATION_NONE; if (startupParam.apotization.BeginsWith("weak", TString::kIgnoreCase)) apodTag = F_APODIZATION_WEAK; else if (startupParam.apotization.BeginsWith("medium", TString::kIgnoreCase)) @@ -919,6 +1131,8 @@ int main(int argc, char *argv[]) else if (startupParam.apotization.BeginsWith("strong", TString::kIgnoreCase)) apodTag = F_APODIZATION_STRONG; + cout << endl << "debug> apodTag = " << apodTag << endl; + double start = millitime(); for (unsigned int i=0; i<fourier.size(); i++) { fourier[i]->Transform(apodTag); From 2d2594c903ef619bd3dcd7499fac43c170ed58e0 Mon Sep 17 00:00:00 2001 From: Suter Andreas <andreas.suter@psi.ch> Date: Fri, 13 Feb 2015 16:56:19 +0100 Subject: [PATCH 3/5] first almost feature complete version of musrFT --- ChangeLog | 8 +- src/classes/Makefile.am | 5 + src/classes/PFourierCanvas.cpp | 1312 +++++++++++++++++++++++++++ src/classes/PMusrCanvas.cpp | 203 +---- src/classes/PPrepFourier.cpp | 53 +- src/include/PFourier.h | 2 + src/include/PFourierCanvas.h | 161 ++++ src/include/PFourierCanvasLinkDef.h | 39 + src/include/PMusrCanvas.h | 4 +- src/include/PPrepFourier.h | 13 +- src/musrFT.cpp | 370 ++++++-- src/musrview.cpp | 5 +- 12 files changed, 1900 insertions(+), 275 deletions(-) create mode 100644 src/classes/PFourierCanvas.cpp create mode 100644 src/include/PFourierCanvas.h create mode 100644 src/include/PFourierCanvasLinkDef.h diff --git a/ChangeLog b/ChangeLog index 4f92c0c7..fcceb54f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,10 +4,10 @@ changes since 0.13.0 =================================== -NEW 2015-02-02 first partial implementation of a standalone Fourier transform/plotter: - musrFT. A LOT of the functionality is still missing! Initially it is - meant to be used for HAL-9500, i.e. Fourier transform WITHOUT lifetime - correction. +NEW 2015-02-13 first implementation of a standalone Fourier transform/plotter: + musrFT. Initially it is meant to be used for HAL-9500, + i.e. Fourier transform WITHOUT lifetime correction. + A first simple minded lifetime correction is implemented as well. NEW 2014-12-18 first implementation of a GLOBAL block which allows to shorten a typical msr-file. Duplicate entries from the RUN blocks can be added here. Furthermore, the 'lifetimecorrection' flag is diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am index 0be0160e..16ede5e0 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -4,6 +4,7 @@ h_sources = \ ../include/PFitterFcn.h \ ../include/PFitter.h \ ../include/PFourier.h \ + ../include/PFourierCanvas.h \ ../include/PFunctionGrammar.h \ ../include/PFunction.h \ ../include/PFunctionHandler.h \ @@ -27,6 +28,7 @@ h_sources_userFcn = \ ../include/PUserFcnBase.h h_linkdef = \ + ../include/PFourierCanvasLinkDef.h \ ../include/PMusrCanvasLinkDef.h \ ../include/PMusrT0LinkDef.h \ ../include/PStartupHandlerLinkDef.h @@ -35,6 +37,7 @@ h_linkdef_userFcn = \ ../include/PUserFcnBaseLinkDef.h dict_h_sources = \ + PFourierCanvasDict.h \ PMusrCanvasDict.h \ PMusrT0Dict.h \ PStartupHandlerDict.h @@ -46,6 +49,7 @@ cpp_sources = \ PFitter.cpp \ PFitterFcn.cpp \ PFourier.cpp \ + PFourierCanvas.cpp \ PFunction.cpp \ PFunctionHandler.cpp \ PMsr2Data.cpp \ @@ -68,6 +72,7 @@ cpp_sources_userFcn = \ PUserFcnBase.cpp dict_cpp_sources = \ + PFourierCanvasDict.cpp \ PMusrCanvasDict.cpp \ PMusrT0Dict.cpp \ PStartupHandlerDict.cpp diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp new file mode 100644 index 00000000..a8466779 --- /dev/null +++ b/src/classes/PFourierCanvas.cpp @@ -0,0 +1,1312 @@ +/*************************************************************************** + + PFourierCanvas.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2015 by Andreas Suter * + * andreas.suter@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#include <iostream> +#include <fstream> +using namespace std; + +#include <TColor.h> +#include <TRandom.h> +#include <TROOT.h> +#include <TDatime.h> +#include <TMath.h> + +#include "PFourierCanvas.h" + +#define YINFO 0.2 +#define YTITLE 0.95 + +ClassImpQ(PFourierCanvas) + +//--------------------------------------------------------------------------- +/** + * <p>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; +} + +//--------------------------------------------------------------------------- +/** + * <p>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<PFourier*> &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 +} + +//--------------------------------------------------------------------------- +/** + * <p>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<PFourier*> &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 +} + +//--------------------------------------------------------------------------- +/** + * <p>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; i<fFourierHistos.size(); i++) { + delete fFourierHistos[i].dataFourierRe; + delete fFourierHistos[i].dataFourierIm; + delete fFourierHistos[i].dataFourierPwr; + delete fFourierHistos[i].dataFourierPhase; + } + fFourierHistos.clear(); + } + + CleanupAverage(); + +/* + if (fFourierPad) { + fFourierPad->Clear(); + delete fFourierPad; + fFourierPad = 0; + } +*/ + if (fInfoPad) { + fInfoPad->Clear(); + delete fInfoPad; + fInfoPad = 0; + } + +/* + if (fMainCanvas) { + delete fMainCanvas; + fMainCanvas = 0; + } +*/ +} + +//-------------------------------------------------------------------------- +// Done (SIGNAL) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>Filters keyboard events, and if they are a command key (see below) carries out the + * necessary actions. + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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; i<fFourier.size(); i++) { + strncpy(title, fFourier[i]->GetDataTitle(), sizeof(title)); + // add entry + fInfoPad->AddEntry(fFourierHistos[i].dataFourierRe, title, "p"); + } + + fInfoPad->Draw(); + fMainCanvas->cd(); + fMainCanvas->Update(); +} + +//-------------------------------------------------------------------------- +// SetTimeout (public) +//-------------------------------------------------------------------------- +/** + * <p>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 +//-------------------------------------------------------------------------- +/** + * <p>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 +//-------------------------------------------------------------------------- +/** + * <p>Saves the currently displayed Fourier data set in ascii column format. + */ +void PFourierCanvas::SaveDataAscii() +{ + // STILL MISSING +} + +//-------------------------------------------------------------------------- +// CreateXaxisTitle (private) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p> 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) +//-------------------------------------------------------------------------- +/** + * <p>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; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierRe = fFourier[i]->GetRealFourier(); + 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; i<fFourierHistos.size(); i++) { + for (int j=start; j<=end; j++) { + dval = fFourierHistos[i].dataFourierRe->GetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=1; j<fFourierHistos[i].dataFourierRe->GetNbinsX(); j++) { + fFourierHistos[i].dataFourierRe->SetBinContent(j, fFourierHistos[i].dataFourierRe->GetBinContent(j)/fabs(max)); + } + } + // imaginary + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=start; j<=end; j++) { + dval = fFourierHistos[i].dataFourierIm->GetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=1; j<fFourierHistos[i].dataFourierIm->GetNbinsX(); j++) { + fFourierHistos[i].dataFourierIm->SetBinContent(j, fFourierHistos[i].dataFourierIm->GetBinContent(j)/fabs(max)); + } + } + // power + max = 0.0; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=start; j<=end; j++) { + dval = fFourierHistos[i].dataFourierPwr->GetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=1; j<fFourierHistos[i].dataFourierPwr->GetNbinsX(); j++) { + fFourierHistos[i].dataFourierPwr->SetBinContent(j, fFourierHistos[i].dataFourierPwr->GetBinContent(j)/fabs(max)); + } + } + // phase + max = 0.0; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=start; j<=end; j++) { + dval = fFourierHistos[i].dataFourierPhase->GetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + } + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + for (int j=1; j<fFourierHistos[i].dataFourierPhase->GetNbinsX(); j++) { + fFourierHistos[i].dataFourierPhase->SetBinContent(j, fFourierHistos[i].dataFourierPhase->GetBinContent(j)/fabs(max)); + } + } + + // set the marker and line color + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + if (i < fColorList.size()) { + fFourierHistos[i].dataFourierRe->SetMarkerColor(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; i<fFourierHistos.size(); i++) { + if (i < fMarkerList.size()) { + fFourierHistos[i].dataFourierRe->SetMarkerStyle(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) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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; j<fFourierHistos[0].dataFourierRe->GetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + if (j < fFourierHistos[i].dataFourierRe->GetNbinsX()) + 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; j<fFourierHistos[0].dataFourierIm->GetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + if (j < fFourierHistos[i].dataFourierIm->GetNbinsX()) + 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; j<fFourierHistos[0].dataFourierPwr->GetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + if (j < fFourierHistos[i].dataFourierPwr->GetNbinsX()) + 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; j<fFourierHistos[0].dataFourierPhase->GetNbinsX(); j++) { + dval = 0.0; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + if (j < fFourierHistos[i].dataFourierPhase->GetNbinsX()) + 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) +//-------------------------------------------------------------------------- +/** + * <p>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<fFourierHistos.size(); i++) { + if (GetMaximum(fFourierHistos[i].dataFourierRe, xmin, xmax) > 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; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierRe->Draw("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<fFourierHistos.size(); i++) { + if (GetMaximum(fFourierHistos[i].dataFourierIm, xmin, xmax) > 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; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierIm->Draw("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<fFourierHistos.size(); i++) { + if (GetMaximum(fFourierHistos[i].dataFourierRe, xmin, xmax) > 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; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierRe->Draw("psame"); + } + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierIm->Draw("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<fFourierHistos.size(); i++) { + if (GetMaximum(fFourierHistos[i].dataFourierPwr, xmin, xmax) > 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; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierPwr->Draw("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<fFourierHistos.size(); i++) { + if (GetMaximum(fFourierHistos[i].dataFourierPhase, xmin, xmax) > 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; i<fFourierHistos.size(); i++) { + fFourierHistos[i].dataFourierPhase->Draw("psame"); + } + break; + default: + break; + } + + fFourierPad->Update(); + + fMainCanvas->cd(); + fMainCanvas->Update(); +} + +//-------------------------------------------------------------------------- +// PlotFourierPhaseValue (private) +//-------------------------------------------------------------------------- +/** + * <p>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) +//-------------------------------------------------------------------------- +/** + * <p>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("<Real>"); + 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("<Imag>"); + 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("<Real+Imag>"); + 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("<Power>"); + fFourierAverage.dataFourierPwr->GetXaxis()->SetTitle(fXaxisTitle.Data()); + fFourierAverage.dataFourierPwr->Draw("p"); + break; + case FOURIER_PLOT_PHASE: + fFourierAverage.dataFourierPhase->GetYaxis()->SetTitle("<Phase>"); + 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) +//-------------------------------------------------------------------------- +/** + * <p>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; i<fFourierHistos.size(); i++) { // loop over all data sets + if ((fFourierHistos[i].dataFourierRe != 0) && (fFourierHistos[i].dataFourierIm != 0)) { + for (Int_t j=0; j<fFourierHistos[i].dataFourierRe->GetNbinsX(); 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) +//-------------------------------------------------------------------------- +/** + * <p>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; i<fFourierHistos.size(); i++) { // loop over all data sets + if ((fFourierHistos[i].dataFourierRe != 0) && (fFourierHistos[i].dataFourierIm != 0)) { + for (Int_t j=0; j<fFourierHistos[i].dataFourierRe->GetNbinsX(); 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) +//-------------------------------------------------------------------------- +/** + * <p>returns the maximum of a histogram in the range [xmin, xmax]. + * If xmin = xmax = -1.0, the whole histogram range is used. + * + * <b>return:</b> + * - 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; i<end; i++) { + binContent = histo->GetBinContent(i); + if (max < binContent) + max = binContent; + } + + return max; +} + +//-------------------------------------------------------------------------- +// GetMinimum (private) +//-------------------------------------------------------------------------- +/** + * <p>returns the minimum of a histogram in the range [xmin, xmax]. + * If xmin = xmax = -1.0, the whole histogram range is used. + * + * <b>return:</b> + * - 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; i<end; i++) { + binContent = histo->GetBinContent(i); + if (min > binContent) + min = binContent; + } + + return min; +} + +//-------------------------------------------------------------------------- +// GetInterpolatedValue (private) +//-------------------------------------------------------------------------- +/** + * <p>search for xVal in histo. If xVal is not found exactly, interpolate and + * return the interpolated y-value. + * + * <b>return:</b> + * - interpolated value if xVal is within histo range, 0 otherwise. + * + * \param histo pointer of the histogram + * \param xVal x-value to be looked for + */ +Double_t PFourierCanvas::GetInterpolatedValue(TH1F* histo, Double_t xVal) +{ + if (histo == 0) + return 0.0; + + Int_t idx = histo->FindBin(xVal); + + // make sure idx is within range + if ((idx < 1) || (idx > histo->GetNbinsX())) + return 0.0; + + // make sure the lower bound idx is found. This is needed since + // FindBin rounds towards the closer idx. + if (histo->GetBinCenter(idx) > xVal) + --idx; + + Double_t x0, x1, y0, y1; + x0 = histo->GetBinCenter(idx); + x1 = histo->GetBinCenter(idx+1); + y0 = histo->GetBinContent(idx); + y1 = histo->GetBinContent(idx+1); + + return (y1-y0)*(xVal-x0)/(x1-x0)+y0; +} diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 9762e17b..4fcf2bd5 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -3797,15 +3797,10 @@ void PMusrCanvas::HandleAverage() // calculate all the average data sets double dval; if (fDataAvg.data != 0) { - if (!CalcAlignment(eTime)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: data: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].data->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].data->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].data, fData[0].data->GetBinContent(i)); } fDataAvg.data->SetBinContent(i, dval/fData.size()); } @@ -3816,15 +3811,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.data->SetMarkerStyle(fData[0].data->GetMarkerStyle()); } if (fDataAvg.dataFourierRe != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Re: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].dataFourierRe->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].dataFourierRe->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierRe, fData[0].dataFourierRe->GetBinContent(i)); } fDataAvg.dataFourierRe->SetBinContent(i, dval/fData.size()); } @@ -3835,15 +3825,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierRe->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle()); } if (fDataAvg.dataFourierIm != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Im: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].dataFourierIm->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].dataFourierIm->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierIm, fData[0].dataFourierIm->GetBinContent(i)); } fDataAvg.dataFourierIm->SetBinContent(i, dval/fData.size()); } @@ -3854,15 +3839,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierIm->SetMarkerStyle(fData[0].dataFourierIm->GetMarkerStyle()); } if (fDataAvg.dataFourierPwr != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Pwr: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].dataFourierPwr->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].dataFourierPwr->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierPwr, fData[0].dataFourierPwr->GetBinContent(i)); } fDataAvg.dataFourierPwr->SetBinContent(i, dval/fData.size()); } @@ -3873,15 +3853,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierPwr->SetMarkerStyle(fData[0].dataFourierPwr->GetMarkerStyle()); } if (fDataAvg.dataFourierPhase != 0) { - if (!CalcAlignment(eFreq)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Phase: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].dataFourierPhase->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].dataFourierPhase->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].dataFourierPhase, fData[0].dataFourierPhase->GetBinContent(i)); } fDataAvg.dataFourierPhase->SetBinContent(i, dval/fData.size()); } @@ -3892,15 +3867,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); } if (fDataAvg.theory != 0) { - if (!CalcAlignment(eTheoTime)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: theory: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].theory->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].theory->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].theory, fData[0].theory->GetBinContent(i)); } fDataAvg.theory->SetBinContent(i, dval/fData.size()); } @@ -3910,7 +3880,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].theoryFourierRe->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].theoryFourierRe->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierRe, fData[0].theoryFourierRe->GetBinContent(i)); } fDataAvg.theoryFourierRe->SetBinContent(i, dval/fData.size()); } @@ -3924,7 +3894,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].theoryFourierIm->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].theoryFourierIm->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierIm, fData[0].theoryFourierIm->GetBinContent(i)); } fDataAvg.theoryFourierIm->SetBinContent(i, dval/fData.size()); } @@ -3938,7 +3908,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].theoryFourierPwr->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].theoryFourierPwr->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierPwr, fData[0].theoryFourierPwr->GetBinContent(i)); } fDataAvg.theoryFourierPwr->SetBinContent(i, dval/fData.size()); } @@ -3952,7 +3922,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].theoryFourierPhase->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].theoryFourierPhase->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].theoryFourierPhase, fData[0].theoryFourierPhase->GetBinContent(i)); } fDataAvg.theoryFourierPhase->SetBinContent(i, dval/fData.size()); } @@ -3963,15 +3933,10 @@ void PMusrCanvas::HandleAverage() fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle()); } if (fDataAvg.diff != 0) { - if (!CalcAlignment(eTime)) { - cerr << endl << ">> PMusrCanvas::HandleAverage: diff: **WARNING** only approx. alignment possible." << endl; - } for (Int_t i=0; i<fData[0].diff->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - if ((i-fAlignmentOffset[j]) > 0) { - dval += fData[j].diff->GetBinContent(i-fAlignmentOffset[j]); - } + dval += GetInterpolatedValue(fData[j].diff, fData[0].diff->GetBinContent(i)); } fDataAvg.diff->SetBinContent(i, dval/fData.size()); } @@ -3985,7 +3950,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].diffFourierRe->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].diffFourierRe->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierRe, fData[0].diffFourierRe->GetBinContent(i)); } fDataAvg.diffFourierRe->SetBinContent(i, dval/fData.size()); } @@ -3999,7 +3964,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].diffFourierIm->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].diffFourierIm->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierIm, fData[0].diffFourierIm->GetBinContent(i)); } fDataAvg.diffFourierIm->SetBinContent(i, dval/fData.size()); } @@ -4013,7 +3978,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].diffFourierPwr->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].diffFourierPwr->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierPwr, fData[0].diffFourierPwr->GetBinContent(i)); } fDataAvg.diffFourierPwr->SetBinContent(i, dval/fData.size()); } @@ -4027,7 +3992,7 @@ void PMusrCanvas::HandleAverage() for (Int_t i=0; i<fData[0].diffFourierPhase->GetNbinsX(); i++) { dval = 0.0; for (UInt_t j=0; j<fData.size(); j++) { - dval += fData[j].diffFourierPhase->GetBinContent(i); + dval += GetInterpolatedValue(fData[j].diffFourierPhase, fData[0].diffFourierPhase->GetBinContent(i)); } fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size()); } @@ -6314,130 +6279,40 @@ UInt_t PMusrCanvas::GetNeededAccuracy(PMsrParamStructure param) //-------------------------------------------------------------------------- -// CalcAlignment (private) +// GetInterpolatedValue (private) //-------------------------------------------------------------------------- /** - * <p>Calculates the alignment index for each data set needed to average the data. + * <p>search for xVal in histo. If xVal is not found exactly, interpolate and + * return the interpolated y-value. * * <b>return:</b> - * - 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; i<fData.size(); i++) { - if (fData[i].data->GetXaxis()->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; i<fData.size(); i++) { - dval = fData[i].data->GetXaxis()->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 && (++j<fData[idx].data->GetNbinsX())); - } - } 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; i<fData.size(); i++) { - if (fData[i].theory->GetXaxis()->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; i<fData.size(); i++) { - dval = fData[i].theory->GetXaxis()->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 && (++j<fData[idx].theory->GetNbinsX())); - } - } 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; i<fData.size(); i++) { - if (fData[i].dataFourierRe->GetXaxis()->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; i<fData.size(); i++) { - dval = fData[i].dataFourierRe->GetXaxis()->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 && (++j<fData[idx].dataFourierRe->GetNbinsX())); - } - } 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; i<fData.size(); i++) { - if (fData[i].theoryFourierRe->GetXaxis()->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; i<fData.size(); i++) { - dval = fData[i].theoryFourierRe->GetXaxis()->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 && (++j<fData[idx].theoryFourierRe->GetNbinsX())); - } - } - - return result; + return (y1-y0)*(xVal-x0)/(x1-x0)+y0; } + diff --git a/src/classes/PPrepFourier.cpp b/src/classes/PPrepFourier.cpp index 30750f8f..cf077229 100644 --- a/src/classes/PPrepFourier.cpp +++ b/src/classes/PPrepFourier.cpp @@ -70,14 +70,12 @@ PPrepFourier::~PPrepFourier() // SetBkgRange //-------------------------------------------------------------------------- /** - * <p> + * <p>set the background range. * - * \param bkgRange + * \param bkgRange array with background range */ void PPrepFourier::SetBkgRange(const Int_t *bkgRange) { - cout << endl << "debug> bkgRange: " << bkgRange[0] << ", " << bkgRange[1] << endl; - int err=0; if (bkgRange[0] >= -1) { fBkgRange[0] = bkgRange[0]; @@ -118,9 +116,9 @@ void PPrepFourier::SetBkgRange(const Int_t *bkgRange) // SetPacking //-------------------------------------------------------------------------- /** - * <p> + * <p>set the packing for the histograms. * - * \param packing + * \param packing number to be used. */ void PPrepFourier::SetPacking(const Int_t packing) { @@ -135,9 +133,10 @@ void PPrepFourier::SetPacking(const Int_t packing) // AddData //-------------------------------------------------------------------------- /** - * <p> + * <p>add a data-set (time domain data + meta information) to the internal + * data vector. * - * \param data + * \param data set to be added */ void PPrepFourier::AddData(musrFT_data &data) { @@ -148,13 +147,10 @@ void PPrepFourier::AddData(musrFT_data &data) // DoBkgCorrection //-------------------------------------------------------------------------- /** - * <p> - * + * <p>Correct the internal data sets according to a background interval given. */ void PPrepFourier::DoBkgCorrection() { - cout << endl << "debug> DoBkgCorrection ..."; - // make sure fData are already present, and if not create the necessary data sets if (fData.size() != fRawData.size()) { InitData(); @@ -162,7 +158,6 @@ void PPrepFourier::DoBkgCorrection() // if no bkg-range is given, nothing needs to be done if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) { - cout << endl << "debug> no background range given ..."; return; } @@ -181,7 +176,6 @@ void PPrepFourier::DoBkgCorrection() bkg += fRawData[i].rawData[j]; } bkg /= (fBkgRange[1]-fBkgRange[0]+1); - cout << endl << "debug> histo: " << i << ", bkg=" << bkg; // correct data for (unsigned int j=0; j<fData[i].size(); j++) @@ -193,13 +187,10 @@ void PPrepFourier::DoBkgCorrection() // DoPacking //-------------------------------------------------------------------------- /** - * <p> - * + * <p>Rebin (pack) the internal data. */ void PPrepFourier::DoPacking() { - cout << endl << "debug> DoPacking : pack=" << fPacking << " ..."; - // make sure fData are already present, and if not create the necessary data sets if (fData.size() != fRawData.size()) { InitData(); @@ -224,10 +215,6 @@ void PPrepFourier::DoPacking() // change the original data set with the packed one fData[i].clear(); fData[i] = tmpData; - - cout << endl << "debug> histo " << i+1 << ": packed data: "; - for (unsigned int j=0; j<15; j++) - cout << fData[i][j] << ", "; } } @@ -235,27 +222,25 @@ void PPrepFourier::DoPacking() // DoFiltering //-------------------------------------------------------------------------- /** - * <p> - * + * <p>Not implemented yet. */ void PPrepFourier::DoFiltering() { - cout << endl << "debug> DoFiltering not yet implemented ..."; - // make sure fData are already present, and if not create the necessary data sets if (fData.size() != fRawData.size()) { InitData(); } - } //-------------------------------------------------------------------------- // DoLifeTimeCorrection //-------------------------------------------------------------------------- /** - * <p> + * <p>Try to do a muon life time correction. The idea is to estimate N0 without + * any theory. This will be OK for high fields (> couple kGauss) but not so good + * for low fields. * - * \param fudge + * \param fudge rescaling factor for the estimated N0. Should be around 1 */ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge) { @@ -294,9 +279,9 @@ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge) // GetInfo //-------------------------------------------------------------------------- /** - * <p> + * <p>Returns the meta information of a data set. * - * \param idx + * \param idx index of the object */ TString PPrepFourier::GetInfo(const UInt_t idx) { @@ -312,8 +297,8 @@ TString PPrepFourier::GetInfo(const UInt_t idx) // GetData //-------------------------------------------------------------------------- /** - * <p> - * + * <p>Creates the requested TH1F objects and returns them. The ownership is with + * the caller. */ vector<TH1F*> PPrepFourier::GetData() { @@ -382,7 +367,7 @@ vector<TH1F*> PPrepFourier::GetData() * <p>Creates the requested TH1F object and returns it. The ownership is with * the caller. * - * \param idx + * \param idx index of the requested histogram */ TH1F *PPrepFourier::GetData(const UInt_t idx) { diff --git a/src/include/PFourier.h b/src/include/PFourier.h index acb7a5b2..787450e1 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -52,6 +52,8 @@ class PFourier virtual void Transform(UInt_t apodizationTag = 0); + virtual const char* GetDataTitle() { return fData->GetTitle(); } + virtual const Int_t GetUnitTag() { return fUnitTag; } virtual Double_t GetResolution() { return fResolution; } virtual TH1F* GetRealFourier(const Double_t scale = 1.0); virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0); diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h new file mode 100644 index 00000000..0634fb96 --- /dev/null +++ b/src/include/PFourierCanvas.h @@ -0,0 +1,161 @@ +/*************************************************************************** + * Copyright (C) 2007-2015 by Andreas Suter * + * andreas.suter@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#ifndef _PFOURIERCANVAS_H_ +#define _PFOURIERCANVAS_H_ + +#include <TObject.h> +#include <TQObject.h> +#include <TTimer.h> +#include <TStyle.h> +#include <TRootCanvas.h> +#include <TGMenu.h> +#include <TCanvas.h> +#include <TPaveText.h> +#include <TLegend.h> +#include <TPad.h> +#include <TH1F.h> +#include <TLatex.h> + +#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 + +//------------------------------------------------------------------------ +/** + * <p>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; + +//------------------------------------------------------------------------ +/** + * <p>typedef to make to code more readable: list of histogram data sets. + */ +typedef vector<PFourierCanvasDataSet> PFourierCanvasDataList; + +//-------------------------------------------------------------------------- +/** + * <p> + */ +class PFourierCanvas : public TObject, public TQObject +{ + public: + PFourierCanvas(); + PFourierCanvas(vector<PFourier*> &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<PFourier*> &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<PFourier*> fFourier; ///< keeps all the Fourier data, ownership is with the caller + PFourierCanvasDataList fFourierHistos; ///< keeps all the Fourier histos + PFourierCanvasDataSet fFourierAverage; ///< keeps the average of the Fourier histos + Double_t fCurrentFourierPhase; ///< keeps the current Fourier phase (real/imag) + TLatex *fCurrentFourierPhaseText; ///< used in Re/Im Fourier to show the current phase in the pad + + TStyle *fStyle; ///< A collection of all graphics attributes + + TTimer *fTimeoutTimer; ///< timeout timer in order to terminate if no action is taking place for too long + + PIntVector fMarkerList; ///< list of markers + PIntVector fColorList; ///< list of colors + + // canvas menu related variables + TRootCanvas *fImp; ///< ROOT native GUI version of main window with menubar and drawing area + TGMenuBar *fBar; ///< menu bar + TGPopupMenu *fPopupMain; ///< popup menu MusrFT in the main menu bar + TGPopupMenu *fPopupSave; ///< popup menu of the MusrFT/Save Data sub menu + TGPopupMenu *fPopupFourier; ///< popup menu of the MusrFT/Fourier sub menu + + // canvas related variables + TCanvas *fMainCanvas; ///< main canvas + TPaveText *fTitlePad; ///< title pad used to display a title + TPad *fFourierPad; ///< fourier pad used to display the fourier + TLegend *fInfoPad; ///< info pad used to display a legend of the data plotted + + virtual void CreateXaxisTitle(); + virtual void CreateStyle(); + virtual void InitFourierDataSets(); + virtual void InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh); + virtual void CleanupAverage(); + virtual void HandleAverage(); + + virtual void PlotFourier(); + virtual void PlotFourierPhaseValue(); + virtual void PlotAverage(); + virtual void IncrementFourierPhase(); + virtual void DecrementFourierPhase(); + + virtual Double_t GetMaximum(TH1F* histo, Double_t xmin=-1.0, Double_t xmax=-1.0); + virtual Double_t GetMinimum(TH1F* histo, Double_t xmin=-1.0, Double_t xmax=-1.0); + virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal); + + ClassDef(PFourierCanvas, 1) +}; + +#endif // _PFOURIERCANVAS_H_ diff --git a/src/include/PFourierCanvasLinkDef.h b/src/include/PFourierCanvasLinkDef.h new file mode 100644 index 00000000..8ceb2d47 --- /dev/null +++ b/src/include/PFourierCanvasLinkDef.h @@ -0,0 +1,39 @@ +/*************************************************************************** + + PMusrCanvasLinkDef.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2015 by Andreas Suter * + * andreas.suter@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class PFourierCanvas+; + +#endif + diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 7b44260e..e8296e4c 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -278,9 +278,7 @@ class PMusrCanvas : public TObject, public TQObject PRunListCollection *fRunList; ///< data handler #endif // __MAKECINT__ - enum EAlignTag {eTime, eTheoTime, eFreq, eTheoFreq}; PMusrCanvasDataSet fDataAvg; ///< set of all averaged data to be plotted (asymmetry/single histogram) - PIntVector fAlignmentOffset; ///< holds the vector with the time/freq alignment offsets PMusrCanvasDataList fData; ///< list of all histogram data to be plotted (asymmetry/single histogram) PMusrCanvasNonMusrDataList fNonMusrData; ///< list of all error graphs to be plotted (non-muSR) @@ -333,7 +331,7 @@ class PMusrCanvas : public TObject, public TQObject virtual Bool_t IsScaleN0AndBkg(); virtual UInt_t GetNeededAccuracy(PMsrParamStructure param); - virtual Bool_t CalcAlignment(const EAlignTag tag); + virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal); ClassDef(PMusrCanvas, 1) }; diff --git a/src/include/PPrepFourier.h b/src/include/PPrepFourier.h index 735a23b7..d1066809 100644 --- a/src/include/PPrepFourier.h +++ b/src/include/PPrepFourier.h @@ -40,14 +40,15 @@ using namespace std; //---------------------------------------------------------------------------- /** - * <p> + * <p>Data structure holding raw time domain uSR data together with some + * necessary meta information. */ typedef struct { - TString info; //< keeps all the meta information - double timeResolution; //< time resolution in (usec) - int t0; //< keep the t0 bin - Double_t timeRange[2]; //< time range to be used, given in (usec). - PDoubleVector rawData; //< a single time domain data vector + TString info; ///< keeps all the meta information + double timeResolution; ///< time resolution in (usec) + int t0; ///< keep the t0 bin + Double_t timeRange[2]; ///< time range to be used, given in (usec). + PDoubleVector rawData; ///< a single time domain data vector } musrFT_data; //---------------------------------------------------------------------------- diff --git a/src/musrFT.cpp b/src/musrFT.cpp index 3551eedb..edf84510 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -38,6 +38,8 @@ #include <vector> using namespace std; +#include <TApplication.h> +#include <TROOT.h> #include <TString.h> #include <TObjArray.h> #include <TObjString.h> @@ -50,36 +52,39 @@ using namespace std; #include "PRunDataHandler.h" #include "PPrepFourier.h" #include "PFourier.h" +#include "PFourierCanvas.h" //---------------------------------------------------------------------------- /** - * <p> + * <p>Structure keeping the command line options. */ typedef struct { - vector<TString> msrFln; - vector<TString> dataFln; - vector<TString> dataFileFormat; - TString graphicFormat; - TString dumpFln; - int bkg[2]; - TString fourierOpt; - TString apotization; - int fourierPower; - TString fourierUnits; - double initialPhase; - double fourierRange[2]; - double timeRange[2]; - vector<int> histo; - bool showAverage; - vector<int> t0; - int packing; - TString title; - double lifetimecorrection; //< is == 0.0 for NO life time correction, otherwise it holds the fudge factor + vector<TString> msrFln; ///< msr-file names to be used. + vector<TString> dataFln; ///< raw-data-file names to be used. + vector<TString> 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<int> 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<int> 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; //------------------------------------------------------------------------- /** - * <p> + * <p>prints the musrFT usage. */ void musrFT_syntax() { @@ -104,9 +109,9 @@ void musrFT_syntax() cout << endl << " --filter : filter and filter-specific-information -- ***TO BE WRITTEN YET***."; cout << endl << " -b, --background <start> <end>: background interval used to estimate the backround to be"; cout << endl << " subtracted before the Fourier transform. <start>, <end> to be given in bins."; - cout << endl << " -fo, --fourier-option <fopt>: <fopt> can be 'real', 'power', 'imag', 'real+imag', of 'phase'."; + cout << endl << " -fo, --fourier-option <fopt>: <fopt> can be 'real', 'imag', 'real+imag', 'power', or 'phase'."; cout << endl << " If this is not defined (neither on the command line nor in the musrFT_startup.xml),"; - cout << endl << " default will be 'real'."; + cout << endl << " default will be 'power'."; cout << endl << " -apod, --apodization <val> : <val> can be either 'none', 'weak', 'medium', 'strong'."; cout << endl << " Default will be 'none'."; cout << endl << " -fp, --fourier-power <N> : <N> being the Fourier power, i.e. 2^<N> used for zero padding."; @@ -115,7 +120,7 @@ void musrFT_syntax() cout << endl << " One may choose between the fields (Gauss) or (Tesla), the frequency (MHz),"; cout << endl << " and the angular-frequency domain (Mc/s)."; cout << endl << " Default will be 'MHz'."; - cout << endl << " -ph, --phase <val> : defines the inital phase <val>. This only is of concern for 'real',"; + cout << endl << " -ph, --phase <val> : defines the initial phase <val>. This only is of concern for 'real',"; cout << endl << " '<imag>', and 'real+imag'."; cout << endl << " Default will be 0.0."; cout << endl << " -fr, --fourier-range <start> <end> : Fourier range. <start>, <end> are interpreted in the units given."; @@ -134,12 +139,14 @@ void musrFT_syntax() cout << endl << " to all histos."; cout << endl << " Example: musrFT -df lem15_his_01234.root -fo real --t0 2750 --histo 1 3"; cout << endl << " -pa, --packing <N> : if <N> (an integer), the time domain data will first be packed/rebinned by <N>."; - cout << endl << " -t, --title <title> : give a global title for the plot."; + cout << endl << " --title <title> : give a global title for the plot."; cout << endl << " --create-msr-file <fln> : creates a msr-file based on the command line options"; cout << endl << " provided. This will help on the way to a full fitting model."; cout << endl << " ***TO BE WRITTEN YET.***"; cout << endl << " -lc, --lifetimecorrection <fudge>: try to eliminate muon life time decay. Only makes sense for low"; cout << endl << " transverse fields. <fudge> is a tweaking factor and should be kept around 1.0."; + cout << endl << " --timeout <timeout> : <timeout> given in seconds after which musrFT terminates."; + cout << endl << " If <timeout> <= 0, no timeout will take place. Default <timeout> is 3600."; cout << endl << endl; } @@ -147,16 +154,17 @@ void musrFT_syntax() /** * <p>initialize startup parameters. * - * \param startupParam + * \param startupParam command line options */ void musrFT_init(musrFT_startup_param &startupParam) { startupParam.graphicFormat = TString(""); startupParam.dumpFln = TString(""); + startupParam.msrFlnOut = TString(""); startupParam.bkg[0] = -1; startupParam.bkg[1] = -1; - startupParam.fourierOpt = TString("real"); - startupParam.apotization = TString("none"); + startupParam.fourierOpt = TString("??"); + startupParam.apodization = TString("none"); startupParam.fourierPower = -1; startupParam.fourierUnits = TString("??"); startupParam.initialPhase = 0.0; @@ -168,6 +176,7 @@ void musrFT_init(musrFT_startup_param &startupParam) startupParam.packing = 1; startupParam.title = TString(""); startupParam.lifetimecorrection = 0.0; + startupParam.timeout = 3600; } //------------------------------------------------------------------------- @@ -273,9 +282,9 @@ bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &s * * <b>return:</b> 0 if everything is OK, 1 for --version or --help, 2 for an error. * - * \param argc - * \param argv - * \param startupParam + * \param argc number of command line arguments + * \param argv command line argument array + * \param startupParam command line data structure */ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupParam) { @@ -356,7 +365,7 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa cerr << endl << ">> musrFT **ERROR** found option --apodization with unrecognized argument '" << topt << "'." << endl; return 2; } - startupParam.apotization = topt; + startupParam.apodization = topt; i++; } else if (tstr.BeginsWith("-fp") || tstr.BeginsWith("--fourier-power")) { if (i+1 >= argc) { // something is wrong since there needs to be two arguments here @@ -449,7 +458,7 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa cerr << endl << ">> musrFT **ERROR** found option --t0 without argument!" << endl; return 2; } - } else if (tstr.BeginsWith("-t ") || tstr.BeginsWith("--title")) { + } else if (tstr.BeginsWith("--title")) { if (i+1 >= argc) { // something is wrong since there needs to be an argument here cerr << endl << ">> musrFT **ERROR** found option --title without argument!" << endl; return 2; @@ -468,6 +477,13 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa return 2; } startupParam.packing = pack.Atoi(); + } else if (tstr.BeginsWith("--create-msr-file")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --create-msr-file without argument!" << endl; + return 2; + } + ++i; + startupParam.msrFlnOut = TString(argv[i]); } else if (tstr.BeginsWith("-lc") || tstr.BeginsWith("--lifetimecorrection")) { if (i+1 >= argc) { // something is wrong since there needs to be an argument here cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection without argument!" << endl; @@ -476,10 +492,22 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa ++i; TString fudge(argv[i]); if (!fudge.IsFloat()) { - cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection with a fudge which is not an double '" << fudge << "'." << endl; + cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection with a fudge which is not a double '" << fudge << "'." << endl; return 2; } startupParam.lifetimecorrection = fudge.Atof(); + } else if (tstr.BeginsWith("--timeout")) { + if (i+1 >= argc) { // something is wrong since there needs to be an argument here + cerr << endl << ">> musrFT **ERROR** found option --timeout without argument!" << endl; + return 2; + } + ++i; + TString tt(argv[i]); + if (!tt.IsDigit()) { + cerr << endl << ">> musrFT **ERROR** found option --timeout with a <timeout> which is not an integer '" << tt << "'." << endl; + return 2; + } + startupParam.timeout = tt.Atoi(); } else if (tstr.BeginsWith("-df") || tstr.BeginsWith("--data-file")) { while (++i < argc) { if (argv[i][0] == '-') { @@ -541,7 +569,11 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa //---------------------------------------------------------------------------------------- /** + * <p>Collects the meta information form the raw-data-file. * + * \param fln file name of the raw-data-file + * \param rawRunData raw-data-file object + * \param metaInfo return string which will contain the meta information. */ void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &metaInfo) { @@ -581,6 +613,13 @@ void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &met } //------------------------------------------------------------------------- +/** + * <p>Estimates the t0's of the raw-data-files. It simply is looking for the + * maximum of the raw-data (assuming a prompt peak). This will fail for LEM + * and ISIS data for sure. + * + * \param rd raw-data-file collection (see PPrepFourier.h) + */ void musrFT_estimateT0(musrFT_data &rd) { cout << endl << ">> musrFT **WARNING** try to estimate t0 from maximum in the data set"; @@ -600,12 +639,12 @@ void musrFT_estimateT0(musrFT_data &rd) } //------------------------------------------------------------------------- -void musrFT_cleanup(TH1F *h) /** - * <p> + * <p> deletes a histogram. * - * \param h histogram to be deleted + * \param h point to a ROOT histogram object */ +void musrFT_cleanup(TH1F *h) { if (h) { delete h; @@ -615,9 +654,12 @@ void musrFT_cleanup(TH1F *h) //------------------------------------------------------------------------- /** - * <p> + * <p>Dump the Fourier transformed data into an ascii file. * - * \param data + * \param fln dump file name + * \param fourierData collection of all the Fourier transformed data. + * \param start starting point from where the data shall be written to file. + * \param end ending point up to where the data shall be written to file. */ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, double end) { @@ -691,14 +733,15 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d return 0; } -//--------------------------------------------------- +//------------------------------------------------------------------------- /** - * <p> + * <p>Groups the histograms before Fourier transform. This is used to group + * detectors. * - * \param runDataHandler - * \param global - * \param run - * \param rd + * \param runDataHandler raw-run-data object containing the data + * \param global pointer to the GLOBAL block of the msr-file + * \param run reference to the relevant RUN block of the msr-file + * \param rd data collection which will hold the grouped histograms. */ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, PMsrRunBlock &run, musrFT_data &rd) { @@ -779,7 +822,115 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, return 0; } -//--------------------------------------------------- +//------------------------------------------------------------------------- +/** + * <p>Dumps an msr-file according to the given command line settings. This + * is meant to generate an initial msr-file for a given data-file. This + * routine is 'stupid' in the sense that it knows nothing about the data-files. + * Hence when feeding it with senseless command line settings, the resulting + * msr-file fed back to musrFT might do funny things! + * + * \param param command line options + */ +void musrFT_dumpMsrFile(musrFT_startup_param ¶m) +{ + ofstream fout(param.msrFlnOut.Data(), ofstream::out); + + // write title + if (param.title.Length() == 0) { // create title if not given + if (param.dataFln.size() != 0) { + param.title = param.dataFln[0]; + } else { + param.title = param.msrFlnOut; + } + } + fout << param.title << endl; + fout << "###############################################################" << endl; + + // write GLOBAL block + fout << "GLOBAL" << endl; + fout << "fittype 0 (single histogram fit)" << endl; + if (param.t0.size() == 1) { // only a single t0 value given, hence assume it is valid for ALL histos + fout << "t0 " << param.t0[0] << endl; + } + if ((param.timeRange[0] != -1.0) && (param.timeRange[1] != -1.0)) { + fout << "fit " << param.timeRange[0] << " " << param.timeRange[1] << endl; + } + fout << "packing " << param.packing << endl; + fout << endl; + fout << "###############################################################" << endl; + + // write RUN block + // get extension of the data file + TString fileFormat("MUSR-ROOT"); + for (unsigned int i=0; i<param.dataFln.size(); i++) { + if (param.dataFileFormat[i].BeginsWith("PsiBin")) + fileFormat = TString("PSI-MDU"); + else if (param.dataFileFormat[i].BeginsWith("NeXus")) + fileFormat = TString("NEXUS"); + else if (param.dataFileFormat[i].BeginsWith("Mud")) + fileFormat = TString("MUD"); + for (unsigned int j=0; j<param.histo.size(); j++) { + fout << "RUN " << param.dataFln[i] << " BXXX IXX " << fileFormat << " (name beamline institute data-file-format)" << endl; + fout << "forward " << param.histo[j] << endl; + if ((param.t0.size() > 1) && (j < param.t0.size())) { + fout << "t0 " << param.t0[j] << endl; + } + if ((param.bkg[0] > -1) && (param.bkg[1] > -1)) + fout << "background " << param.bkg[0] << " " << param.bkg[1] << endl; + fout << "#--------------------------------------------------------------" << endl; + } + } + fout << endl; + fout << "###############################################################" << endl; + + // write PLOT block + fout << "PLOT 0 (single histo plot)" << endl; + if (param.histo.size() == 0) { + fout << "runs 1" << endl; + } else { + fout << "runs "; + for (unsigned int i=0; i<param.histo.size(); i++) + fout << i+1 << " "; + fout << endl; + } + if ((param.timeRange[0] == -1.0) && (param.timeRange[1] == -1.0)) { + fout << "range 0 10" << endl; + } else { + fout << "range " << param.timeRange[0] << " " << param.timeRange[1] << endl; + } + fout << endl; + fout << "###############################################################" << endl; + + // write FOURIER block + fout << "FOURIER" << endl; + if (param.fourierUnits.BeginsWith("??")) { // Fourier units not given, hence choose MHz + fout << "units MHz # units either 'Gauss', 'MHz', or 'Mc/s'" << endl; + } else { + fout << "units " << param.fourierUnits << " # units either 'Gauss', 'MHz', or 'Mc/s'" << endl; + } + if (param.fourierOpt.BeginsWith("??")) { // Fourier plot option not given, hence choose POWER + fout << "plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE" << endl; + } else { + fout << "plot " << param.fourierOpt << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE" << endl; + } + if (param.fourierPower > 1) { + fout << "fourier_power " << param.fourierPower << endl; + } + fout << "apodization " << param.apodization << " # NONE, WEAK, MEDIUM, STRONG" << endl; + if ((param.fourierRange[0] > -1.0) && (param.fourierRange[1] > -1.0)) { + fout << "range " << param.fourierRange[0] << " " << param.fourierRange[1] << endl; + } + + fout.close(); +} + +//------------------------------------------------------------------------- +/** + * <p>Gets time a time stamp in msec. Used to measure the calculation time. + * + * <b>return:</b> time stamp with msec resolution. + */ double millitime() { struct timeval now; @@ -790,15 +941,20 @@ double millitime() //------------------------------------------------------------------------- /** - * <p> + * <p>musrFT is used to do a Fourier transform of uSR data without any fitting. + * It directly Fourier transforms the raw histogram data (exception see --lifetimecorrection), + * and hence only will give staisfactory results for applied fields of larger a + * couple of kGauss. It is meant to be used to get a feeling what time-domain + * model will be appropriate. It is NOT meant for ANY quantitative analysis! * - * \param argc - * \param argv + * \param argc number of command line arguments + * \param argv command line argument array */ int main(int argc, char *argv[]) { Int_t unitTag = FOURIER_UNIT_NOT_GIVEN; Int_t apodTag = F_APODIZATION_NONE; + Int_t fourierPlotTag = FOURIER_PLOT_NOT_GIVEN; // only program name alone if (argc == 1) { @@ -821,6 +977,12 @@ int main(int argc, char *argv[]) return retVal; } + // dump msr-file + if (startupParam.msrFlnOut.Length() > 0) { + musrFT_dumpMsrFile(startupParam); + return PMUSR_SUCCESS; + } + // read startup file char startup_path_name[128]; PStartupOptions startup_options; @@ -903,7 +1065,6 @@ int main(int argc, char *argv[]) // load data-file(s) provided directly for (unsigned int i=msrHandler.size(); i<msrHandler.size()+startupParam.dataFln.size(); i++) { -// cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()]; // create run data handler if (startupHandler) runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList()); @@ -956,6 +1117,9 @@ int main(int argc, char *argv[]) unsigned int idx=0; // get meta info, time resolution, time range, raw data sets if (i < msrHandler.size()) { // obtain info from msr-files + // keep title if not overwritten by the command line + if (startupParam.title.Length() == 0) + startupParam.title = *(msrHandler[0]->GetMsrTitle()); // keep PLOT block info PMsrPlotList *plot = msrHandler[i]->GetMsrPlotList(); if (plot == 0) { @@ -993,7 +1157,11 @@ int main(int argc, char *argv[]) // get range if ((startupParam.fourierRange[0] == -1) && (startupParam.fourierRange[1] == -1)) { // no Fourier range given from the command line startupParam.fourierRange[0] = fourierBlock->fPlotRange[0]; - startupParam.fourierRange[1] = fourierBlock->fPlotRange[1]; + startupParam.fourierRange[1] = fourierBlock->fPlotRange[1]; + } + // get Fourier plot option, i.e. real, imag, power, phase + if (startupParam.fourierOpt.BeginsWith("??")) { // only do something if not overwritten by the command line + fourierPlotTag = fourierBlock->fPlotTag; } } } @@ -1030,8 +1198,10 @@ int main(int argc, char *argv[]) rd.timeRange[0] = startupParam.timeRange[0]; rd.timeRange[1] = startupParam.timeRange[1]; } else { - rd.timeRange[0] = plot->at(0).fTmin[0]; - rd.timeRange[1] = plot->at(0).fTmax[0]; + if (plot->at(0).fTmin.size() > 0) { + rd.timeRange[0] = plot->at(0).fTmin[0]; + rd.timeRange[1] = plot->at(0).fTmax[0]; + } } // handle data set(s) @@ -1089,6 +1259,22 @@ int main(int argc, char *argv[]) } } + // make sure Fourier plot tag is set + if (fourierPlotTag == FOURIER_PLOT_NOT_GIVEN) { + if (!startupParam.fourierOpt.CompareTo("real", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_REAL; + else if (!startupParam.fourierOpt.CompareTo("imag", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_IMAG; + else if (!startupParam.fourierOpt.CompareTo("real+imag", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_REAL_AND_IMAG; + else if (!startupParam.fourierOpt.CompareTo("power", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_POWER; + else if (!startupParam.fourierOpt.CompareTo("phase", TString::kIgnoreCase)) + fourierPlotTag = FOURIER_PLOT_PHASE; + else + fourierPlotTag = FOURIER_PLOT_POWER; + } + // calculate background levels and subtract them from the data data.DoBkgCorrection(); @@ -1124,15 +1310,13 @@ int main(int argc, char *argv[]) } // Fourier transform data - if (startupParam.apotization.BeginsWith("weak", TString::kIgnoreCase)) + if (startupParam.apodization.BeginsWith("weak", TString::kIgnoreCase)) apodTag = F_APODIZATION_WEAK; - else if (startupParam.apotization.BeginsWith("medium", TString::kIgnoreCase)) + else if (startupParam.apodization.BeginsWith("medium", TString::kIgnoreCase)) apodTag = F_APODIZATION_MEDIUM; - else if (startupParam.apotization.BeginsWith("strong", TString::kIgnoreCase)) + else if (startupParam.apodization.BeginsWith("strong", TString::kIgnoreCase)) apodTag = F_APODIZATION_STRONG; - cout << endl << "debug> apodTag = " << apodTag << endl; - double start = millitime(); for (unsigned int i=0; i<fourier.size(); i++) { fourier[i]->Transform(apodTag); @@ -1140,19 +1324,83 @@ int main(int argc, char *argv[]) double end = millitime(); cout << endl << "debug> after FFT. calculation time: " << (end-start)/1.0e3 << " (sec)." << endl; + PFourierCanvas *fourierCanvas = 0; + // if Fourier dumped if whished do it now if (startupParam.dumpFln.Length() > 0) { musrFT_dumpData(startupParam.dumpFln, fourier, startupParam.fourierRange[0], startupParam.fourierRange[1]); - } else { + } else { // do Canvas - // if Fourier graphical export is whished, switch to batch mode + // if Fourier graphical export is whished, switch to batch mode + Bool_t batch = false; + if (startupParam.graphicFormat.Length() != 0) { + batch = true; + argv[argc] = (char*)malloc(16*sizeof(char)); + strcpy(argv[argc], "-b"); + argc++; + } - // plot the Fourier transform + // plot the Fourier transform + TApplication app("App", &argc, argv); + if (startupHandler) { + fourierCanvas = new PFourierCanvas(fourier, startupParam.title.Data(), + startupParam.showAverage, fourierPlotTag, + startupParam.fourierRange, startupParam.initialPhase, + 10, 10, 800, 800, + startupHandler->GetMarkerList(), + startupHandler->GetColorList(), + batch); + } else { + fourierCanvas = new PFourierCanvas(fourier, startupParam.title.Data(), + startupParam.showAverage, fourierPlotTag, + startupParam.fourierRange, startupParam.initialPhase, + 10, 10, 800, 800, + batch); + } + + fourierCanvas->UpdateFourierPad(); + fourierCanvas->UpdateInfoPad(); + + Bool_t ok = true; + if (!fourierCanvas->IsValid()) { + cerr << endl << ">> musrFT **SEVERE ERROR** Couldn't invoke all necessary objects, will quit."; + cerr << endl; + ok = false; + } else { + // connect signal/slot + TQObject::Connect("TCanvas", "Closed()", "PFourierCanvas", fourierCanvas, "LastCanvasClosed()"); + + fourierCanvas->SetTimeout(startupParam.timeout); + + fourierCanvas->Connect("Done(Int_t)", "TApplication", &app, "Terminate(Int_t)"); + + if (startupParam.graphicFormat.Length() != 0) { + TString fileName(""); + // create output filename based on the msr- or raw-data-filename + if (startupParam.dataFln.size() > 0) { + fileName = startupParam.dataFln[0]; + } + if (startupParam.msrFln.size() > 0) { + fileName = startupParam.msrFln[0]; + } + Ssiz_t idx = fileName.Last('.'); + fileName.Remove(idx, fileName.Length()); + fileName += "."; + fileName += startupParam.graphicFormat; + fourierCanvas->SaveGraphicsAndQuit(fileName.Data()); + } + } + // check that everything is ok + if (ok) + app.Run(true); // true needed that Run will return after quit so that cleanup works } // cleanup + if (fourierCanvas) + delete fourierCanvas; + if (startupHandler) delete startupHandler; diff --git a/src/musrview.cpp b/src/musrview.cpp index 6e035945..17d3d978 100644 --- a/src/musrview.cpp +++ b/src/musrview.cpp @@ -340,9 +340,9 @@ int main(int argc, char *argv[]) } if (asciiOutput) { - // save data in batch mode + // save data in batch mode musrCanvas->SaveDataAscii(); - musrCanvas->Done(0); + musrCanvas->Done(0); } // keep musrCanvas objects @@ -360,7 +360,6 @@ int main(int argc, char *argv[]) sprintf(canvasName, "fMainCanvas%d", i); if (gROOT->GetListOfCanvases()->FindObject(canvasName) != 0) { canvasVector[i]->~PMusrCanvas(); - } else { } } canvasVector.empty(); From f9197cb815f8f1380d959fed91a3012e546eee4d Mon Sep 17 00:00:00 2001 From: Suter Andreas <andreas.suter@psi.ch> Date: Fri, 13 Feb 2015 17:04:50 +0100 Subject: [PATCH 4/5] removed outdated help information. --- src/musrFT.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/musrFT.cpp b/src/musrFT.cpp index edf84510..68d3c47b 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -142,7 +142,6 @@ void musrFT_syntax() cout << endl << " --title <title> : give a global title for the plot."; cout << endl << " --create-msr-file <fln> : creates a msr-file based on the command line options"; cout << endl << " provided. This will help on the way to a full fitting model."; - cout << endl << " ***TO BE WRITTEN YET.***"; cout << endl << " -lc, --lifetimecorrection <fudge>: try to eliminate muon life time decay. Only makes sense for low"; cout << endl << " transverse fields. <fudge> is a tweaking factor and should be kept around 1.0."; cout << endl << " --timeout <timeout> : <timeout> given in seconds after which musrFT terminates."; From e097cc3119e5e98400a6949dbabb58661facb2c4 Mon Sep 17 00:00:00 2001 From: Andreas Suter <andreas.suter@psi.ch> Date: Sat, 14 Feb 2015 17:46:06 +0100 Subject: [PATCH 5/5] added data export from the GUI menu. --- src/classes/PFourierCanvas.cpp | 196 ++++++++++++++++++++++++++++++--- src/include/PFourierCanvas.h | 8 +- 2 files changed, 185 insertions(+), 19 deletions(-) diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp index a8466779..9c33bd04 100644 --- a/src/classes/PFourierCanvas.cpp +++ b/src/classes/PFourierCanvas.cpp @@ -36,12 +36,17 @@ using namespace std; #include <TROOT.h> #include <TDatime.h> #include <TMath.h> +#include <TGFileDialog.h> #include "PFourierCanvas.h" #define YINFO 0.2 #define YTITLE 0.95 +static const char *gFiletypes[] = { "All files", "*", + "Data files", "*.dat", + 0, 0 }; + ClassImpQ(PFourierCanvas) //--------------------------------------------------------------------------- @@ -70,7 +75,6 @@ PFourierCanvas::PFourierCanvas() fImp = 0; fBar = 0; fPopupMain = 0; - fPopupSave = 0; fPopupFourier = 0; fMainCanvas = 0; @@ -170,8 +174,6 @@ PFourierCanvas::PFourierCanvas(vector<PFourier*> &fourier, const Char_t* title, */ PFourierCanvas::~PFourierCanvas() { - cout << endl << "debug> in ~PFourierCanvas()." << endl; - if (fTimeoutTimer) delete fTimeoutTimer; @@ -345,8 +347,16 @@ void PFourierCanvas::HandleMenuPopup(Int_t id) CleanupAverage(); PlotFourier(); } - } else if (id == P_MENU_ID_SAVE_DATA+P_MENU_ID_SAVE_ASCII) { - cout << endl << "debug> Save still missing ..." << endl; + } else if (id == P_MENU_ID_EXPORT_DATA) { + static TString dir("."); + TGFileInfo fi; + fi.fFileTypes = gFiletypes; + fi.fIniDir = StrDup(dir); + fi.fOverwrite = true; + new TGFileDialog(0, fImp, kFDSave, &fi); + if (fi.fFilename && strlen(fi.fFilename)) { + ExportData(fi.fFilename); + } } } @@ -460,14 +470,176 @@ void PFourierCanvas::SaveGraphicsAndQuit(const Char_t *fileName) } //-------------------------------------------------------------------------- -// SaveDataAscii +// ExportData //-------------------------------------------------------------------------- /** - * <p>Saves the currently displayed Fourier data set in ascii column format. + * <p>Exports the currently displayed Fourier data set in ascii column format. */ -void PFourierCanvas::SaveDataAscii() +void PFourierCanvas::ExportData(const Char_t *pathFileName) { - // STILL MISSING + TString pfn(""); + + if (pathFileName) { // path file name provided + pfn = TString(pathFileName); + } else { // path file name NOT provided, generate a default path file name + cerr << endl << ">> PFourierCanvas::ExportData **ERROR** NO path file name provided. Will do nothing." << endl; + return; + } + + TString xAxis(""), yAxis(""); + Int_t xMinBin, xMaxBin; + if (fAveragedView) { + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + xAxis = fFourierHistos[0].dataFourierRe->GetXaxis()->GetTitle(); + yAxis = TString("<Real>"); + xMinBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetLast(); + break; + case FOURIER_PLOT_IMAG: + xAxis = fFourierHistos[0].dataFourierIm->GetXaxis()->GetTitle(); + yAxis = TString("<Imag>"); + xMinBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetLast(); + break; + case FOURIER_PLOT_POWER: + xAxis = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetTitle(); + yAxis = TString("<Power>"); + xMinBin = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetLast(); + break; + case FOURIER_PLOT_PHASE: + xAxis = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetTitle(); + yAxis = TString("<Phase>"); + xMinBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetLast(); + break; + default: + xAxis = TString("??"); + yAxis = TString("??"); + xMinBin = 0; + xMaxBin = 0; + break; + } + } else { + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + xAxis = fFourierHistos[0].dataFourierRe->GetXaxis()->GetTitle(); + yAxis = TString("Real"); + xMinBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetLast(); + break; + case FOURIER_PLOT_IMAG: + xAxis = fFourierHistos[0].dataFourierIm->GetXaxis()->GetTitle(); + yAxis = TString("Imag"); + xMinBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetLast(); + break; + case FOURIER_PLOT_POWER: + xAxis = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetTitle(); + yAxis = TString("Power"); + xMinBin = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetLast(); + break; + case FOURIER_PLOT_PHASE: + xAxis = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetTitle(); + yAxis = TString("Phase"); + xMinBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetLast(); + break; + default: + xAxis = TString("??"); + yAxis = TString("??"); + xMinBin = 0; + xMaxBin = 0; + break; + } + } + + // write data to file + ofstream fout(pfn.Data(), ofstream::out); + + if (fAveragedView) { + // write header + fout << "% " << pfn << endl; + fout << "% averaged data of:" << endl; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + fout << "% " << fFourierHistos[i].dataFourierRe->GetTitle() << endl; + } + fout << "%------------" << endl; + fout << "% " << xAxis << ", " << yAxis << endl; + for (int i=xMinBin; i<xMaxBin; i++) { + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + fout << fFourierAverage.dataFourierRe->GetBinCenter(i) << ", " << fFourierAverage.dataFourierRe->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_IMAG: + fout << fFourierAverage.dataFourierIm->GetBinCenter(i) << ", " << fFourierAverage.dataFourierIm->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_POWER: + fout << fFourierAverage.dataFourierPwr->GetBinCenter(i) << ", " << fFourierAverage.dataFourierPwr->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_PHASE: + fout << fFourierAverage.dataFourierPhase->GetBinCenter(i) << ", " << fFourierAverage.dataFourierPhase->GetBinContent(i) << endl; + break; + default: + break; + } + } + } else { + // write header + fout << "% " << pfn << endl; + fout << "% data of:" << endl; + for (unsigned int i=0; i<fFourierHistos.size(); i++) { + fout << "% " << fFourierHistos[i].dataFourierRe->GetTitle() << endl; + } + fout << "%------------" << endl; + fout << "% "; + for (unsigned int i=0; i<fFourierHistos.size()-1; i++) { + fout << xAxis << i << ", " << yAxis << i << ", "; + } + fout << xAxis << fFourierHistos.size()-1 << ", " << yAxis << fFourierHistos.size()-1 << endl; + + // write data + for (int i=xMinBin; i<xMaxBin; i++) { + for (unsigned int j=0; j<fFourierHistos.size()-1; j++) { + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + fout << fFourierHistos[j].dataFourierRe->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierRe->GetBinContent(i) << ", "; + break; + case FOURIER_PLOT_IMAG: + fout << fFourierHistos[j].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierIm->GetBinContent(i) << ", "; + break; + case FOURIER_PLOT_POWER: + fout << fFourierHistos[j].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPwr->GetBinContent(i) << ", "; + break; + case FOURIER_PLOT_PHASE: + fout << fFourierHistos[j].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhase->GetBinContent(i) << ", "; + break; + default: + break; + } + } + switch (fCurrentPlotView) { + case FOURIER_PLOT_REAL: + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_IMAG: + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_POWER: + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_PHASE: + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinContent(i) << endl; + break; + default: + break; + } + } + } + + fout.close(); } //-------------------------------------------------------------------------- @@ -664,7 +836,6 @@ void PFourierCanvas::InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t w fImp = 0; fBar = 0; fPopupMain = 0; - fPopupSave = 0; fPopupFourier = 0; fMainCanvas = 0; @@ -726,10 +897,7 @@ void PFourierCanvas::InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t w 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); + fPopupMain->AddEntry("Export Data", P_MENU_ID_EXPORT_DATA); fBar->MapSubwindows(); fBar->Layout(); diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h index 0634fb96..d509e39e 100644 --- a/src/include/PFourierCanvas.h +++ b/src/include/PFourierCanvas.h @@ -40,7 +40,7 @@ // 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_EXPORT_DATA 10003 #define P_MENU_ID_FOURIER_REAL 100 #define P_MENU_ID_FOURIER_IMAG 101 @@ -50,8 +50,6 @@ #define P_MENU_ID_FOURIER_PHASE_PLUS 105 #define P_MENU_ID_FOURIER_PHASE_MINUS 106 -#define P_MENU_ID_SAVE_ASCII 200 - //------------------------------------------------------------------------ /** * <p>Structure holding all necessary Fourier histograms. @@ -99,7 +97,7 @@ class PFourierCanvas : public TObject, public TQObject virtual void SetTimeout(Int_t ival); virtual void SaveGraphicsAndQuit(const Char_t *fileName); - virtual void SaveDataAscii(); + virtual void ExportData(const Char_t *pathFileName); private: Int_t fTimeout; ///< timeout after which the Done signal should be emited. If timeout <= 0, no timeout is taking place @@ -129,7 +127,7 @@ class PFourierCanvas : public TObject, public TQObject 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 *fPopupSave; ///< popup menu of the MusrFT/Save Data sub menu TGPopupMenu *fPopupFourier; ///< popup menu of the MusrFT/Fourier sub menu // canvas related variables