From f7834aaead4696e2e295c7fb821bff66e8b75e14 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Tue, 29 Nov 2016 11:16:45 +0100 Subject: [PATCH] added class PFTPhaseCorrection which allows to calculate the phase optimzied real part of the Fourier transform. musrFT is already using it, for musrview it is still missing. --- src/classes/PFourier.cpp | 287 +++++++++++++++++++++++++++++---- src/classes/PFourierCanvas.cpp | 55 ++++--- src/include/PFourier.h | 55 ++++++- src/include/PFourierCanvas.h | 2 +- 4 files changed, 339 insertions(+), 60 deletions(-) diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index ec516096..b6034b26 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -37,7 +37,9 @@ using namespace std; #include "TF1.h" #include "TAxis.h" -//#include "TFile.h" +#include "Minuit2/FunctionMinimum.h" +#include "Minuit2/MnUserParameters.h" +#include "Minuit2/MnMinimize.h" #include "PMusr.h" #include "PFourier.h" @@ -45,6 +47,227 @@ using namespace std; #define PI 3.14159265358979312 #define PI_HALF 1.57079632679489656 +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// PFTPhaseCorrection +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + * + */ +PFTPhaseCorrection::PFTPhaseCorrection(const Int_t minBin, const Int_t maxBin) : + fMinBin(minBin), fMaxBin(maxBin) +{ + Init(); +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + * + */ +PFTPhaseCorrection::PFTPhaseCorrection(vector &reFT, vector &imFT, const Int_t minBin, const Int_t maxBin) : + fReal(reFT), fImag(imFT), fMinBin(minBin), fMaxBin(maxBin) +{ + Init(); + + if (fMinBin == -1) + fMinBin = 0; + if (fMaxBin == -1) + fMaxBin = fReal.size(); + if (fMaxBin > fReal.size()) + fMaxBin = fReal.size(); + + fRealPh.resize(fReal.size()); + fImagPh.resize(fReal.size()); +} + +//-------------------------------------------------------------------------- +// Minimize (public) +//-------------------------------------------------------------------------- +/** + * + */ +void PFTPhaseCorrection::Minimize() +{ + // create Minuit2 parameters + ROOT::Minuit2::MnUserParameters upar; + + upar.Add("c0", fPh_c0, 0.5); + upar.Add("c1", fPh_c1, 0.5); + + // create minimizer + ROOT::Minuit2::MnMinimize mn_min(*this, upar); + + // minimize + ROOT::Minuit2::FunctionMinimum min = mn_min(); + + if (min.IsValid()) { + fPh_c0 = min.UserState().Value("c0"); + fPh_c1 = min.UserState().Value("c1"); + fMin = min.Fval(); + } else { + fMin = -1.0; + fValid = false; + cout << endl << ">> **WARNING** minimize failed to find a minimum for the real FT phase correction ..." << endl; + return; + } +} + +//-------------------------------------------------------------------------- +// GetPhaseCorrectionParam (public) +//-------------------------------------------------------------------------- +/** + * + */ +Double_t PFTPhaseCorrection::GetPhaseCorrectionParam(UInt_t idx) +{ + Double_t result=0.0; + + if (idx == 0) + result = fPh_c0; + else if (idx == 1) + result = fPh_c1; + else + cerr << ">> **ERROR** requested phase correction parameter with index=" << idx << " does not exist!" << endl; + + return result; +} + +//-------------------------------------------------------------------------- +// GetMinimum (public) +//-------------------------------------------------------------------------- +/** + * + */ +Double_t PFTPhaseCorrection::GetMinimum() +{ + if (!fValid) { + cerr << ">> **ERROR** requested minimum is invalid!" << endl; + return -1.0; + } + + return fMin; +} + +//-------------------------------------------------------------------------- +// Init (private) +//-------------------------------------------------------------------------- +/** + * + */ +void PFTPhaseCorrection::Init() +{ + fValid = true; + fPh_c0 = 0.0; + fPh_c1 = 0.0; + fGamma = 1.0; + fMin = -1.0; +} + +//-------------------------------------------------------------------------- +// CalcPhasedFT (private) +//-------------------------------------------------------------------------- +/** + * + */ +void PFTPhaseCorrection::CalcPhasedFT() const +{ + Double_t phi=0.0; + Double_t w=0.0; + + for (UInt_t i=0; i 1.0e-15) { + hh = dval / norm; + entropy -= hh * log(hh); + } + } + + return entropy; +} + +//-------------------------------------------------------------------------- +// operator() (private) +//-------------------------------------------------------------------------- +/** + * + */ +double PFTPhaseCorrection::operator()(const vector &par) const +{ + // par : [0]: c0, [1]: c1 + + fPh_c0 = par[0]; + fPh_c1 = par[1]; + + CalcPhasedFT(); + CalcRealPhFTDerivative(); + + return Entropy()+Penalty(); +} + +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// PFourier //-------------------------------------------------------------------------- // Constructor //-------------------------------------------------------------------------- @@ -73,6 +296,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi fValid = true; fIn = 0; fOut = 0; + fPhCorrectedReFT = 0; fApodization = F_APODIZATION_NONE; @@ -160,6 +384,8 @@ PFourier::~PFourier() fftw_free(fIn); if (fOut) fftw_free(fOut); + if (fPhCorrectedReFT) + delete fPhCorrectedReFT; } //-------------------------------------------------------------------------- @@ -262,26 +488,19 @@ TH1F* PFourier::GetRealFourier(const Double_t scale) // GetPhaseOptRealFourier (public) //-------------------------------------------------------------------------- /** - *

returns the phase corrected real Fourier transform. The correction angle is - * past back as well. - *

Currently it simply does the following thing: (i) rotate the complex Fourier - * transform through all angles in 1/2° steps, i.e. - * \f$ f_{\rm rot} = (f_{\rm real} + i f_{\rm imag}) \exp(- \alpha)\f$, - * hence \f$ f_{\rm rot} = f_{\rm real} \cos(\alpha) + f_{\rm imag} \sin(\alpha)\f$. - * (ii) search the maximum of \f$ sum_{\alpha} {\cal R}\{f_{\rm rot}\}\f$ for all - * \f$\alpha\f$. From this one gets \f$\alpha_{\rm opt}\f$. (iii) The 'optimal' - * real Fourier transform is \f$f_{\rm opt} = (f_{\rm real} + i f_{\rm imag}) - * \exp(- \alpha_{\rm opt})\f$. + *

returns the phase corrected real Fourier transform. * * \return the TH1F histo of the phase 'optimzed' real Fourier transform. * - * \param phase return value of the optimal phase + * \param phase return value of the optimal phase dispersion phase[0]+phase[1]*i/N * \param scale normalisation factor * \param min minimal freq / field from which to optimise. Given in the choosen unit. * \param max maximal freq / field up to which to optimise. Given in the choosen unit. */ -TH1F* PFourier::GetPhaseOptRealFourier(Double_t &phase, const Double_t scale, const Double_t min, const Double_t max) +TH1F* PFourier::GetPhaseOptRealFourier(vector &phase, const Double_t scale, const Double_t min, const Double_t max) { + phase.resize(2); // c0, c1 + UInt_t noOfFourierBins = 0; if (fNoOfBins % 2 == 0) noOfFourierBins = fNoOfBins/2; @@ -312,32 +531,28 @@ TH1F* PFourier::GetPhaseOptRealFourier(Double_t &phase, const Double_t scale, co imagF.push_back(fOut[i][1]); } - // rotate trough real/imag Fourier through 360° with a 1/2° resolution and keep the integral - Double_t rotRealIntegral[720]; - Double_t sum = 0.0; - Double_t da = 8.72664625997164774e-03; // pi / 360 - for (UInt_t i=0; i<720; i++) { - sum = 0.0; - for (UInt_t j=0; jMinimize(); + if (!fPhCorrectedReFT->IsValid()) { + fValid = false; + cerr << endl << "**ERROR** could not fina a valid phase correction minimum." << endl; + return 0; } - - // keep the optimal phase - phase = maxRotBin*da; + phase[0] = fPhCorrectedReFT->GetPhaseCorrectionParam(0); + phase[1] = fPhCorrectedReFT->GetPhaseCorrectionParam(1); // clean up + if (fPhCorrectedReFT) { + delete fPhCorrectedReFT; + fPhCorrectedReFT = 0; + } realF.clear(); imagF.clear(); @@ -355,8 +570,10 @@ TH1F* PFourier::GetPhaseOptRealFourier(Double_t &phase, const Double_t scale, co } // fill realFourier vector + Double_t ph; for (UInt_t i=0; iSetBinContent(i+1, scale*(fOut[i][0]*cos(phase) + fOut[i][1]*sin(phase))); + ph = phase[0] + phase[1] * (Double_t)((Int_t)i-(Int_t)minBin) / (Double_t)(maxBin-minBin); + realPhaseOptFourier->SetBinContent(i+1, scale*(fOut[i][0]*cos(ph) - fOut[i][1]*sin(ph))); realPhaseOptFourier->SetBinError(i+1, 0.0); } diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp index 07813ff4..0b9cbc1f 100644 --- a/src/classes/PFourierCanvas.cpp +++ b/src/classes/PFourierCanvas.cpp @@ -624,6 +624,12 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName) xMinBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetFirst(); xMaxBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetLast(); break; + case FOURIER_PLOT_REAL_AND_IMAG: + xAxis = fFourierHistos[0].dataFourierRe->GetXaxis()->GetTitle(); + yAxis = TString("Real+Imag"); + xMinBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetLast(); + break; case FOURIER_PLOT_POWER: xAxis = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetTitle(); yAxis = TString("Power"); @@ -703,42 +709,46 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName) for (UInt_t j=0; jGetBinCenter(i) << ", " << fFourierHistos[j].dataFourierRe->GetBinContent(i) << ", "; - break; + 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; + fout << fFourierHistos[j].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierIm->GetBinContent(i) << ", "; + break; + case FOURIER_PLOT_REAL_AND_IMAG: + break; case FOURIER_PLOT_POWER: - fout << fFourierHistos[j].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPwr->GetBinContent(i) << ", "; - break; + 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; + fout << fFourierHistos[j].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhase->GetBinContent(i) << ", "; + break; case FOURIER_PLOT_PHASE_OPT_REAL: - fout << fFourierHistos[j].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhaseOptReal->GetBinContent(i) << ", "; - break; + fout << fFourierHistos[j].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhaseOptReal->GetBinContent(i) << ", "; + break; default: - break; + break; } } switch (fCurrentPlotView) { case FOURIER_PLOT_REAL: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinContent(i) << endl; - break; + 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; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_REAL_AND_IMAG: + break; case FOURIER_PLOT_POWER: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinContent(i) << endl; - break; + 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; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinContent(i) << endl; + break; case FOURIER_PLOT_PHASE_OPT_REAL: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinContent(i) << endl; - break; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinContent(i) << endl; + break; default: - break; + break; } } } @@ -805,7 +815,6 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierPwr = fFourier[i]->GetPowerFourier(); fFourierHistos[i].dataFourierPhase = fFourier[i]->GetPhaseFourier(); fFourierHistos[i].dataFourierPhaseOptReal = fFourier[i]->GetPhaseOptRealFourier(fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); - cout << "debug> histo[" << i << "]: opt. phase = " << fFourierHistos[i].optPhase * 180.0 / TMath::Pi() << "°" << endl; } // rescale histo to abs(maximum) == 1 diff --git a/src/include/PFourier.h b/src/include/PFourier.h index bff94450..38f7e7ee 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -30,8 +30,13 @@ #ifndef _PFOURIER_H_ #define _PFOURIER_H_ +#include +using namespace std; + #include "fftw3.h" +#include "Minuit2/FCNBase.h" + #include "PMusr.h" #define F_APODIZATION_NONE 1 @@ -39,6 +44,52 @@ #define F_APODIZATION_MEDIUM 3 #define F_APODIZATION_STRONG 4 +/** + * Re Fourier phase correction + */ +class PFTPhaseCorrection : public ROOT::Minuit2::FCNBase +{ + public: + PFTPhaseCorrection(const Int_t minBin=-1, const Int_t maxBin=-1); + PFTPhaseCorrection(vector &reFT, vector &imFT, const Int_t minBin=-1, const Int_t maxBin=-1); + virtual ~PFTPhaseCorrection() {} + + virtual Bool_t IsValid() { return fValid; } + virtual void Minimize(); + + virtual void SetGamma(const Double_t gamma) { fGamma = gamma; } + virtual void SetPh(const Double_t c0, const Double_t c1) { fPh_c0 = c0; fPh_c1 = c1; CalcPhasedFT(); CalcRealPhFTDerivative(); } + + virtual Double_t GetGamma() { return fGamma; } + virtual Double_t GetPhaseCorrectionParam(UInt_t idx); + virtual Double_t GetMinimum(); + + private: + Bool_t fValid; + + vector fReal; /// original real Fourier data set + vector fImag; /// original imag Fourier data set + mutable vector fRealPh; /// phased real Fourier data set + mutable vector fImagPh; /// phased imag Fourier data set + mutable vector fRealPhD; /// 1st derivative of fRealPh + + Double_t fMinBin; /// minimum bin from Fourier range to be used for the phase correction estimate + Double_t fMaxBin; /// maximum bin from Fourier range to be used for the phase correction estimate + mutable Double_t fPh_c0; /// constant part of the phase dispersion used for the phase correction + mutable Double_t fPh_c1; /// linear part of the phase dispersion used for the phase correction + Double_t fGamma; /// gamma parameter to balance between entropy and penalty + Double_t fMin; /// keeps the minimum of the entropy/penalty minimization + + virtual void Init(); + virtual void CalcPhasedFT() const; + virtual void CalcRealPhFTDerivative() const; + virtual Double_t Penalty() const; + virtual Double_t Entropy() const; + + virtual Double_t Up() const { return 1.0; } + virtual Double_t operator()(const vector&) const; +}; + /** * muSR Fourier class. */ @@ -57,7 +108,7 @@ class PFourier virtual Double_t GetResolution() { return fResolution; } virtual Double_t GetMaxFreq(); virtual TH1F* GetRealFourier(const Double_t scale = 1.0); - virtual TH1F* GetPhaseOptRealFourier(Double_t &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); + virtual TH1F* GetPhaseOptRealFourier(vector &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0); virtual TH1F* GetPowerFourier(const Double_t scale = 1.0); virtual TH1F* GetPhaseFourier(const Double_t scale = 1.0); @@ -85,6 +136,8 @@ class PFourier fftw_complex *fIn; ///< real part of the Fourier transform fftw_complex *fOut; ///< imaginary part of the Fourier transform + PFTPhaseCorrection *fPhCorrectedReFT; + virtual void PrepareFFTwInputData(UInt_t apodizationTag); virtual void ApodizeData(Int_t apodizationTag); }; diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h index 1229b34c..3aff49e9 100644 --- a/src/include/PFourierCanvas.h +++ b/src/include/PFourierCanvas.h @@ -71,7 +71,7 @@ typedef struct { TH1F *dataFourierPwr; ///< power spectrum of the Fourier transform of the data histogram TH1F *dataFourierPhase; ///< phase spectrum of the Fourier transform of the data histogram TH1F *dataFourierPhaseOptReal; ///< phase otpimized real Fourier transform of the data histogram - Double_t optPhase; ///< optimal phase which maximizes the real Fourier + vector optPhase; ///< optimal phase which maximizes the real Fourier } PFourierCanvasDataSet; //------------------------------------------------------------------------