From d29375c9125c8e02e40314903b981de5d95f6858 Mon Sep 17 00:00:00 2001 From: nemu Date: Tue, 16 Dec 2008 08:04:38 +0000 Subject: [PATCH] a small step further towards Fourier transform. PFourier most probably not functional yet. --- src/classes/Makefile.PMusr | 3 +- src/classes/PFourier.cpp | 418 ++++++++++++++++++++++++++++++++++++ src/classes/PMsrHandler.cpp | 8 +- src/classes/PMusrCanvas.cpp | 22 +- src/include/PFourier.h | 97 +++++++++ src/include/PMusr.h | 1 + src/include/PMusrCanvas.h | 12 +- 7 files changed, 552 insertions(+), 9 deletions(-) create mode 100644 src/classes/PFourier.cpp create mode 100644 src/include/PFourier.h diff --git a/src/classes/Makefile.PMusr b/src/classes/Makefile.PMusr index 35e632b0..6e11691f 100644 --- a/src/classes/Makefile.PMusr +++ b/src/classes/Makefile.PMusr @@ -100,6 +100,7 @@ OBJS += PFitterFcn.o OBJS += PFitter.o OBJS += PMusrCanvas.o PMusrCanvasDict.o OBJS += PUserFcnBase.o PUserFcnBaseDict.o +OBJS += PFourier.o EXTOBJS = EXTOBJS += MuSR_td_PSI_bin.o @@ -149,4 +150,4 @@ ifeq ($(OS),LINUX) cp -pv $(PMUSRPATH)/*.h $(ROOTSYS)/include endif -cleaninstall: clean install \ No newline at end of file +cleaninstall: clean install diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp new file mode 100644 index 00000000..71778b4f --- /dev/null +++ b/src/classes/PFourier.cpp @@ -0,0 +1,418 @@ +/*************************************************************************** + + PFourier.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2008 by Andreas Suter * + * andreas.suter@psi.c * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#include + +#include +using namespace std; + +#include "TH1F.h" +#include "TF1.h" + +#include "TFile.h" + +#include "PMusr.h" +#include "PFourier.h" + +#define PI 3.14159265358979312 +#define PI_HALF 1.57079632679489656 + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

+ * + * \dataType tag indicating if data is histogram, asymmetry, ... + * \data vector with the real data + */ +PFourier::PFourier(int dataType, TH1F *data, + double startTime, double endTime, + unsigned int zeroPaddingPower, bool estimateN0AndBkg) : + fDataType(dataType), fData(data), + fStartTime(startTime), fEndTime(endTime), + fZeroPaddingPower(zeroPaddingPower) +{ + // some necessary checks and initialization + if (fData == 0) { + cout << endl << "**ERROR** PFourier::PFourier: no valid data" << endl << endl; + fValid = false; + return; + } + + if ((fStartTime < 0.0) || (fEndTime < 0.0)) { + cout << endl << "**ERROR** PFourier::PFourier: no valid start or end time." << endl << endl; + fValid = false; + return; + } + + fValid = true; + fIn = 0; + fOut = 0; + + fApodization = F_APODIZATION_NONE; + + // calculate time resolution in (ns) + fTimeResolution = fData->GetBinWidth(1) * 1000.0; + + // if endTime == 0 set it to the last time slot + if (fEndTime == 0.0) { + int last = fData->GetNbinsX()-1; + fEndTime = fData->GetBinCenter(last); + } else { + fEndTime *= 1000.0; // us -> ns + } + + fStartTime *= 1000.0; // us -> ns + + // swap start and end time if necessary + if (fStartTime > fEndTime) { + double keep = fStartTime; + fStartTime = fEndTime; + fEndTime = keep; + } + +cout << endl << "dB = " << 1.0/(2.0 * F_GAMMA_BAR_MUON * (fEndTime-fStartTime)) << " (G), Bmax = " << 1.0/(2.0 * F_GAMMA_BAR_MUON * fTimeResolution) << " (G)" << endl; + + // try to estimate N0 and Bkg just out of the raw data + if (estimateN0AndBkg) { + EstimateN0AndBkg(); + } + + // calculate start and end bin + unsigned int start = (unsigned int)(fStartTime/fTimeResolution); + unsigned int end = (unsigned int)(fEndTime/fTimeResolution); + fNoOfData = end-start; + + // check if zero padding is whished + if (fZeroPaddingPower > 0) { + fNoOfBins = static_cast(pow(2.0, static_cast(fZeroPaddingPower))); + } else { + fNoOfBins = fNoOfData; + } + + // allocate necessary memory + fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins); + fOut = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins); + + // check if memory allocation has been successful + if ((fIn == 0) || (fOut == 0)) { + fValid = false; + return; + } + + fFFTwPlan = fftw_plan_dft_1d(fNoOfBins, fIn, fOut, FFTW_FORWARD, FFTW_ESTIMATE); + + if (!fFFTwPlan) { + fValid = false; + } +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

+ * + */ +PFourier::~PFourier() +{ +cout << endl << "in ~PFourier() ..." << endl; + if (fValid) { + fftw_destroy_plan(fFFTwPlan); + fftw_free(fIn); + fftw_free(fOut); + } +} + +//-------------------------------------------------------------------------- +// Transform +//-------------------------------------------------------------------------- +/** + *

+ * + * \param apodizationTag 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod. + */ +void PFourier::Transform(unsigned int apodizationTag) +{ + if (!fValid) + return; + + if (fDataType == F_SINGLE_HISTO) { + PrepareSingleHistoFFTwInputData(apodizationTag); + } else { + PrepareFFTwInputData(apodizationTag); + } + + fftw_execute(fFFTwPlan); +} + +//-------------------------------------------------------------------------- +// GetRealFourier +//-------------------------------------------------------------------------- +/** + *

+ * + * \param realFourier + */ +void PFourier::GetRealFourier(TH1F *realFourier) +{ + // check if valid flag is set + if (!fValid) + return; + + // reallocate realFourier + char name[256]; + char title[256]; + strncpy(name, realFourier->GetName(), sizeof(name)); + strncpy(title, realFourier->GetTitle(), sizeof(title)); + if (realFourier) { + delete realFourier; + realFourier = 0; + } + realFourier = new TH1F(name, title, fNoOfBins, -fFieldResolution/2.0, fNoOfBins*fFieldResolution+fFieldResolution/2.0); + if (realFourier == 0) { + fValid = false; + cout << endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << endl; + return; + } + + // fill realFourier vector + for (unsigned int i=0; iSetBinContent(i+1, fOut[i][0]); + } +} + +//-------------------------------------------------------------------------- +// GetImaginaryFourier +//-------------------------------------------------------------------------- +/** + *

+ * + * \param imaginaryFourier + */ +void PFourier::GetImaginaryFourier(TH1F *imaginaryFourier) +{ + // check if valid flag is set + if (!fValid) + return; + + // reallocate imaginaryFourier + char name[256]; + char title[256]; + strncpy(name, imaginaryFourier->GetName(), sizeof(name)); + strncpy(title, imaginaryFourier->GetTitle(), sizeof(title)); + if (imaginaryFourier) { + delete imaginaryFourier; + imaginaryFourier = 0; + } + imaginaryFourier = new TH1F(name, title, fNoOfBins, -fFieldResolution/2.0, fNoOfBins*fFieldResolution+fFieldResolution/2.0); + if (imaginaryFourier == 0) { + fValid = false; + cout << endl << "**SEVERE ERROR** couldn't allocate memory for the imaginary part of the Fourier transform." << endl; + return; + } + + // fill imaginaryFourier vector + for (unsigned int i=0; iSetBinContent(i+1, fOut[i][1]); + } +} + +//-------------------------------------------------------------------------- +// EstimateN0AndBkg +//-------------------------------------------------------------------------- +/** + *

+ * + */ +void PFourier::EstimateN0AndBkg() +{ + int noOfBins = fData->GetNbinsX(); + + TH1F summHisto("summHisto", "summHisto", noOfBins, + -fTimeResolution/2.0, (noOfBins-1)*fTimeResolution + fTimeResolution/2.0); + + // fill summHisto + double value = 0.0; + for (int i=1; iGetBinContent(i); + summHisto.SetBinContent(i, value); + summHisto.SetBinError(i, sqrt(value)); + } + +cout << endl << ">> max.summHisto=" << summHisto.GetBinContent(noOfBins-1) << endl << endl; + + // define fit function + TF1 *func = new TF1("func", "[0]*(1-TMath::Exp(-x/[1]))+[2]*x", + -fTimeResolution/2.0, (noOfBins-1)*fTimeResolution + + fTimeResolution/2.0); + // parameter 0 ~ N0 tau + func->SetParameter(0, summHisto.GetBinContent(noOfBins-1)); + // parameter 1 == tau + func->FixParameter(1, PMUON_LIFETIME*1000.0); + // parameter 2 ~ + func->SetParameter(2, summHisto.GetBinContent(noOfBins-1)/(PMUON_LIFETIME*1000.0)*0.05); + + // do the fit + summHisto.Fit(func, "0QR"); // 0->no Graph, Q->quite, R->time range from function + + // get out the parameters + double A = func->GetParameter(0); + double B = func->GetParameter(2); + +cout << endl << ">> A=" << A << ", B=" << B; +cout << endl << ">> N0/per bin=" << A/(PMUON_LIFETIME*1000.0)*fTimeResolution << ", per bin=" << B*fTimeResolution << endl << endl; + + fN0 = A/(PMUON_LIFETIME*1000.0)*fTimeResolution; + fBkg = B*fTimeResolution; + + // clean up + if (func) { + delete func; + func = 0; + } +} + +//-------------------------------------------------------------------------- +// PrepareSingleHistoFFTwInputData +//-------------------------------------------------------------------------- +/** + *

+ * + */ +void PFourier::PrepareSingleHistoFFTwInputData(unsigned int apodizationTag) +{ + // 1st fill fIn + unsigned int start = (unsigned int)(fStartTime/fTimeResolution); + for (unsigned int i=0; iGetBinContent(i+start+1); + fIn[i][1] = 0.0; + } + for (unsigned int i=fNoOfData; i + * + */ +void PFourier::PrepareFFTwInputData(unsigned int apodizationTag) +{ + // 1st fill fIn + unsigned int start = (unsigned int)(fStartTime/fTimeResolution); + for (unsigned int i=0; iGetBinContent(i+start+1); + fIn[i][1] = 0.0; + } + for (unsigned int i=fNoOfData; i + * + * \param apodizationTag + */ +void PFourier::ApodizeData(int apodizationTag) { + + const double cweak[3] = { 0.384093, -0.087577, 0.703484 }; + const double cmedium[3] = { 0.152442, -0.136176, 0.983734 }; + const double cstrong[3] = { 0.045335, 0.554883, 0.399782 }; + + double c[5]; + for (unsigned int i=0; i<5; i++) { + c[i] = 0.0; + } + + switch (apodizationTag) { + case F_APODIZATION_NONE: + return; + break; + case F_APODIZATION_WEAK: + c[0] = cweak[0]+cweak[1]+cweak[2]; + c[1] = -(cweak[1]+2.0*cweak[2]); + c[2] = cweak[2]; + break; + case F_APODIZATION_MEDIUM: + c[0] = cmedium[0]+cmedium[1]+cmedium[2]; + c[1] = -(cmedium[1]+2.0*cmedium[2]); + c[2] = cmedium[2]; + break; + case F_APODIZATION_STRONG: + c[0] = cstrong[0]+cstrong[1]+cstrong[2]; + c[1] = -2.0*(cstrong[1]+2.0*cstrong[2]); + c[2] = cstrong[1]+6.0*cstrong[2]; + c[3] = -4.0*cstrong[2]; + c[4] = cstrong[2]; + break; + default: + cout << endl << ">> **ERROR** User Apodization tag " << apodizationTag << " unknown, sorry ..." << endl; + break; + } + + double q; + for (unsigned int i=0; i> in PMsrHandler::HandleFourierEntry, Fourier block present .. fourier.fUnits = FOURIER_UNIT_FIELD; } else if (!str.CompareTo("mhz", TString::kIgnoreCase)) { fourier.fUnits = FOURIER_UNIT_FREQ; + } else if (!str.CompareTo("mc/s", TString::kIgnoreCase)) { + fourier.fUnits = FOURIER_UNIT_CYCLES; } else { error = true; continue; @@ -2008,7 +2012,7 @@ cout << endl << ">> in PMsrHandler::HandleFourierEntry, Fourier block present .. cout << endl << "FOURIER block syntax, parameters in [] are optinal:"; cout << endl; cout << endl << "FOURIER"; - cout << endl << "[units Gauss | MHz]"; + cout << endl << "[units Gauss | MHz | Mc/s]"; cout << endl << "[fourier_power n # n is a number such that zero padding up to 2^n will be used]"; cout << endl << " n=0 means no zero padding"; cout << endl << " 0 <= n <= 20 are allowed values"; diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 8a60649a..7f21c4b3 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -327,7 +327,9 @@ void PMusrCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected) } else if (x == 'd') { HandleDifference(); } else if (x == 'f') { - cout << endl << ">> will show the Fourier transform, to be implemented yet." << endl; + if (fPlotType != MSR_PLOT_NON_MUSR) { + cout << endl << ">> will show the Fourier transform, to be implemented yet." << endl; + } } else if (x == '+') { cout << endl << ">> if Fourier is shown, will add +1° to the Fourier." << endl; } else if (x == '-') { @@ -605,6 +607,14 @@ cout << endl; } // handle data HandleNonMusrDataSet(i, runNo, data); + // disable Fourier menus + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_REAL); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_IMAG); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_REAL_AND_IMAG); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PWR); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS); + fPopupFourier->DisableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS); break; default: fValid = false; @@ -1985,3 +1995,13 @@ void PMusrCanvas::SaveDataDb() cout << endl << ">> Data windows saved in db format ..." << endl; } +//-------------------------------------------------------------------------- +// HandleFourier +//-------------------------------------------------------------------------- +/** + *

+ * + */ +void PMusrCanvas::HandleFourier() +{ +} diff --git a/src/include/PFourier.h b/src/include/PFourier.h new file mode 100644 index 00000000..7832efaa --- /dev/null +++ b/src/include/PFourier.h @@ -0,0 +1,97 @@ +/*************************************************************************** + + PFourier.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + + $Id$ + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2008 by Andreas Suter * + * andreas.suter@psi.c * + * * + * 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 _PFOURIER_H_ +#define _PFOURIER_H_ + +#include "fftw3.h" + +#define F_ESTIMATE_N0_AND_BKG true + +#define F_SINGLE_HISTO_RAW 0 +#define F_SINGLE_HISTO 1 +#define F_ASYMMETRY 2 + +#define F_APODIZATION_NONE 0 +#define F_APODIZATION_WEAK 1 +#define F_APODIZATION_MEDIUM 2 +#define F_APODIZATION_STRONG 3 + +// gamma_muon / (2 pi) = 1.355342e-5 (GHz/G) +#define F_GAMMA_BAR_MUON 1.355342e-5 + +class PFourier +{ + public: + PFourier(int dataType, TH1F *data, + double startTime = 0.0, double endTime = 0.0, + unsigned int zeroPaddingPower = 0, + bool estimateN0AndBkg = false); + virtual ~PFourier(); + + virtual void Transform(unsigned int apodizationTag = 0); + + virtual double GetFieldResolution() { return fFieldResolution; } + virtual void GetRealFourier(TH1F *realFourier); + virtual void GetImaginaryFourier(TH1F *imaginaryFourier); + + virtual bool IsValid() { return fValid; } + + private: + bool fValid; + + int fDataType; ///< 0=Single Histo Raw, 1=Single Histo Life Time Corrected, 2=Asymmetry + int fApodization; ///< 0=none, 1=weak, 2=medium, 3=strong + + double fN0; + double fBkg; + + double fTimeResolution; + double fStartTime; + double fEndTime; + unsigned int fZeroPaddingPower; + double fFieldResolution; + + TH1F *fData; + + unsigned int fNoOfData; + unsigned int fNoOfBins; + fftw_plan fFFTwPlan; + fftw_complex *fIn; + fftw_complex *fOut; + + virtual void PrepareSingleHistoFFTwInputData(unsigned int apodizationTag); + virtual void PrepareFFTwInputData(unsigned int apodizationTag); + virtual void ApodizeData(int apodizationTag); + virtual void EstimateN0AndBkg(); +}; + +#endif // _PFOURIER_H_ diff --git a/src/include/PMusr.h b/src/include/PMusr.h index fdbca2c7..0038fa67 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -98,6 +98,7 @@ using namespace std; #define FOURIER_UNIT_NOT_GIVEN 0 #define FOURIER_UNIT_FIELD 1 #define FOURIER_UNIT_FREQ 2 +#define FOURIER_UNIT_CYCLES 3 #define FOURIER_APOD_NOT_GIVEN 0 #define FOURIER_APOD_NONE 1 diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 520d8103..63fc6e62 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -55,11 +55,12 @@ #define XTHEO 0.75 // Current Plot Views -#define PV_DATA 1 -#define PV_FOURIER_REAL 2 -#define PV_FOURIER_IMAG 3 -#define PV_FOURIER_PWR 4 -#define PV_FOURIER_PHASE 5 +#define PV_DATA 1 +#define PV_FOURIER_REAL 2 +#define PV_FOURIER_IMAG 3 +#define PV_FOURIER_REAL_AND_IMAG 4 +#define PV_FOURIER_PWR 5 +#define PV_FOURIER_PHASE 6 // Canvas menu id's #define P_MENU_ID_FOURIER 10001 @@ -210,6 +211,7 @@ class PMusrCanvas : public TObject, public TQObject virtual void HandleDataSet(unsigned int plotNo, unsigned int runNo, PRunData *data); virtual void HandleNonMusrDataSet(unsigned int plotNo, unsigned int runNo, PRunData *data); virtual void HandleDifference(); + virtual void HandleFourier(); virtual double CalculateDiff(const double x, const double y, TH1F *theo); virtual double CalculateDiff(const double x, const double y, TGraphErrors *theo);