LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 4 Activity

added a first version of a optimized phase correction for the real Fourier transform

Browse Source
This commit is contained in:
suter_a
2016-12-18 10:37:50 +01:00
parent c768c27898
commit 63516fc499
7 changed files with 368 additions and 178 deletions
Download Patch File Download Diff File Expand all files Collapse all files

83
src/classes/PFourier.cpp
View File

@ -94,8 +94,8 @@ void PFTPhaseCorrection::Minimize()
// create Minuit2 parameters
ROOT::Minuit2::MnUserParameters upar;
upar.Add("c0", fPh_c0, 0.5);
upar.Add("c1", fPh_c1, 0.5);
upar.Add("c0", fPh_c0, 2.0);
upar.Add("c1", fPh_c1, 2.0);
// create minimizer
ROOT::Minuit2::MnMinimize mn_min(*this, upar);
@ -296,7 +296,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
fValid = true;
fIn = 0;
fOut = 0;
fPhCorrectedReFT = 0;
//as fPhCorrectedReFT = 0;
fApodization = F_APODIZATION_NONE;
@ -384,8 +384,8 @@ PFourier::~PFourier()
fftw_free(fIn);
if (fOut)
fftw_free(fOut);
if (fPhCorrectedReFT)
delete fPhCorrectedReFT;
//as if (fPhCorrectedReFT)
//as delete fPhCorrectedReFT;
}
//--------------------------------------------------------------------------
@ -485,73 +485,79 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
// GetPhaseOptRealFourier (public)
// GetPhaseOptRealFourier (public, static)
//--------------------------------------------------------------------------
/**
* <p>returns the phase corrected real Fourier transform.
*
* \return the TH1F histo of the phase 'optimzed' real Fourier transform.
*
* \param re real part Fourier histogram
* \param im imaginary part Fourier histogram
* \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(vector<Double_t> &phase, const Double_t scale, const Double_t min, const Double_t max)
TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, vector<Double_t> &phase,
const Double_t scale, const Double_t min, const Double_t max)
{
if ((re == 0) || (im == 0))
return 0;
phase.resize(2); // c0, c1
UInt_t noOfFourierBins = 0;
if (fNoOfBins % 2 == 0)
noOfFourierBins = fNoOfBins/2;
else
noOfFourierBins = (fNoOfBins+1)/2;
TAxis *axis = re->GetXaxis();
UInt_t minBin = 0;
UInt_t maxBin = noOfFourierBins;
Int_t minBin = 1;
Int_t maxBin = axis->GetNbins();
Int_t noOfBins = axis->GetNbins();
Double_t res = axis->GetBinWidth(1);
// check if minimum frequency is given. If yes, get the proper minBin
if (min != -1.0) {
minBin = (UInt_t)(min/fResolution);
minBin = axis->FindFixBin(min);
if ((minBin == 0) || (minBin > maxBin)) {
minBin = 1;
cerr << "**WARNING** minimum frequency/field out of range. Will adopted it." << endl;
}
}
// check if maximum frequency is given. If yes, get the proper maxBin
if (max != -1.0) {
maxBin = (UInt_t)(max/fResolution);
if (maxBin >= noOfFourierBins) {
maxBin = noOfFourierBins;
maxBin = axis->FindFixBin(max);
if ((maxBin == 0) || (maxBin > axis->GetNbins())) {
maxBin = axis->GetNbins();
cerr << "**WARNING** maximum frequency/field out of range. Will adopted it." << endl;
}
}
// copy the real/imag Fourier from min to max
vector<Double_t> realF, imagF;
for (UInt_t i=minBin; i<=maxBin; i++) {
realF.push_back(fOut[i][0]);
imagF.push_back(fOut[i][1]);
for (Int_t i=minBin; i<=maxBin; i++) {
realF.push_back(re->GetBinContent(i));
imagF.push_back(im->GetBinContent(i));
}
// optimize real FT phase
fPhCorrectedReFT = new PFTPhaseCorrection(realF, imagF);
if (fPhCorrectedReFT == 0) {
fValid = false;
PFTPhaseCorrection *phCorrectedReFT = new PFTPhaseCorrection(realF, imagF);
if (phCorrectedReFT == 0) {
cerr << endl << "**SEVERE ERROR** couldn't invoke PFTPhaseCorrection object." << endl;
return 0;
}
fPhCorrectedReFT->Minimize();
if (!fPhCorrectedReFT->IsValid()) {
fValid = false;
cerr << endl << "**ERROR** could not fina a valid phase correction minimum." << endl;
phCorrectedReFT->Minimize();
if (!phCorrectedReFT->IsValid()) {
cerr << endl << "**ERROR** could not find a valid phase correction minimum." << endl;
return 0;
}
phase[0] = fPhCorrectedReFT->GetPhaseCorrectionParam(0);
phase[1] = fPhCorrectedReFT->GetPhaseCorrectionParam(1);
phase[0] = phCorrectedReFT->GetPhaseCorrectionParam(0);
phase[1] = phCorrectedReFT->GetPhaseCorrectionParam(1);
// clean up
if (fPhCorrectedReFT) {
delete fPhCorrectedReFT;
fPhCorrectedReFT = 0;
if (phCorrectedReFT) {
delete phCorrectedReFT;
phCorrectedReFT = 0;
}
realF.clear();
imagF.clear();
@ -559,21 +565,20 @@ TH1F* PFourier::GetPhaseOptRealFourier(vector<Double_t> &phase, const Double_t s
// invoke the real phase optimised histo to be filled. Caller is the owner!
Char_t name[256];
Char_t title[256];
snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", fData->GetName());
snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", fData->GetTitle());
snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", re->GetName());
snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", re->GetTitle());
TH1F *realPhaseOptFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0);
TH1F *realPhaseOptFourier = new TH1F(name, title, noOfBins, -res/2.0, (Double_t)(noOfBins-1)*res+res/2.0);
if (realPhaseOptFourier == 0) {
fValid = false;
cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << endl;
return 0;
}
// fill realFourier vector
Double_t ph;
for (UInt_t i=0; i<noOfFourierBins; i++) {
for (Int_t i=0; i<noOfBins; i++) {
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->SetBinContent(i+1, scale*(re->GetBinContent(i+1)*cos(ph) - im->GetBinContent(i+1)*sin(ph)));
realPhaseOptFourier->SetBinError(i+1, 0.0);
}