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

improved Fourier transform. For details see the ChangeLog

Browse Source
This commit is contained in:
suter_a
2012-11-19 13:27:37 +00:00
parent 11d710bce3
commit 17d7ef4b99
22 changed files with 167 additions and 104 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/classes/PFourier.cpp
+52 -19
View File
@@ -32,6 +32,7 @@
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
#include "TH1F.h"
@@ -57,11 +58,12 @@ using namespace std;
* FOURIER_UNIT_FIELD, FOURIER_UNIT_FREQ, FOURIER_UNIT_CYCLES
* \param startTime start time of the data time window
* \param endTime end time of the data time window
* \param dcCorrected if true, removed DC offset from signal before Fourier transformation, otherwise not
* \param zeroPaddingPower if set to values > 0, there will be zero padding up to 2^zeroPaddingPower
*/
PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTime, UInt_t zeroPaddingPower) :
PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTime, Bool_t dcCorrected, UInt_t zeroPaddingPower) :
fData(data), fUnitTag(unitTag), fStartTime(startTime), fEndTime(endTime),
fZeroPaddingPower(zeroPaddingPower)
fDCCorrected(dcCorrected), fZeroPaddingPower(zeroPaddingPower)
{
// some necessary checks and initialization
if (fData == 0) {
@@ -93,9 +95,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
}
// calculate start and end bin
UInt_t start = (UInt_t)(fStartTime/fTimeResolution);
UInt_t end = (UInt_t)(fEndTime/fTimeResolution);
fNoOfData = end-start;
fNoOfData = (UInt_t)((fEndTime-fStartTime)/fTimeResolution);
// check if zero padding is whished
if (fZeroPaddingPower > 0) {
@@ -112,7 +112,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
Double_t resolution = 1.0/(fTimeResolution*fNoOfBins); // in MHz
switch (fUnitTag) {
case FOURIER_UNIT_FIELD:
fResolution = resolution/F_GAMMA_BAR_MUON;
fResolution = resolution/GAMMA_BAR_MUON;
break;
case FOURIER_UNIT_FREQ:
fResolution = resolution;
@@ -219,7 +219,13 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
snprintf(name, sizeof(name), "%s_Fourier_Re", fData->GetName());
snprintf(title, sizeof(title), "%s_Fourier_Re", fData->GetTitle());
TH1F *realFourier = new TH1F(name, title, fNoOfBins/2, -fResolution/2.0, (Double_t)fNoOfBins/2.0*fResolution+fResolution/2.0);
UInt_t noOfFourierBins = 0;
if (fNoOfBins % 2 == 0)
noOfFourierBins = fNoOfBins/2;
else
noOfFourierBins = (fNoOfBins+1)/2;
TH1F *realFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0);
if (realFourier == 0) {
fValid = false;
cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << endl;
@@ -227,7 +233,7 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
}
// fill realFourier vector
for (UInt_t i=0; i<fNoOfBins/2; i++) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
realFourier->SetBinContent(i+1, scale*fOut[i][0]);
realFourier->SetBinError(i+1, 0.0);
}
@@ -254,7 +260,13 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
snprintf(name, sizeof(name), "%s_Fourier_Im", fData->GetName());
snprintf(title, sizeof(title), "%s_Fourier_Im", fData->GetTitle());
TH1F* imaginaryFourier = new TH1F(name, title, fNoOfBins/2, -fResolution/2.0, (Double_t)fNoOfBins/2.0*fResolution+fResolution/2.0);
UInt_t noOfFourierBins = 0;
if (fNoOfBins % 2 == 0)
noOfFourierBins = fNoOfBins/2;
else
noOfFourierBins = (fNoOfBins+1)/2;
TH1F* imaginaryFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0);
if (imaginaryFourier == 0) {
fValid = false;
cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the imaginary part of the Fourier transform." << endl;
@@ -262,7 +274,7 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
}
// fill imaginaryFourier vector
for (UInt_t i=0; i<fNoOfBins/2; i++) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
imaginaryFourier->SetBinContent(i+1, scale*fOut[i][1]);
imaginaryFourier->SetBinError(i+1, 0.0);
}
@@ -290,7 +302,13 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
snprintf(name, sizeof(name), "%s_Fourier_Pwr", fData->GetName());
snprintf(title, sizeof(title), "%s_Fourier_Pwr", fData->GetTitle());
TH1F* pwrFourier = new TH1F(name, title, fNoOfBins/2, -fResolution/2.0, (Double_t)fNoOfBins/2.0*fResolution+fResolution/2.0);
UInt_t noOfFourierBins = 0;
if (fNoOfBins % 2 == 0)
noOfFourierBins = fNoOfBins/2;
else
noOfFourierBins = (fNoOfBins+1)/2;
TH1F* pwrFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0);
if (pwrFourier == 0) {
fValid = false;
cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the power part of the Fourier transform." << endl;
@@ -298,7 +316,7 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
}
// fill powerFourier vector
for (UInt_t i=0; i<fNoOfBins/2; i++) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
pwrFourier->SetBinContent(i+1, scale*sqrt(fOut[i][0]*fOut[i][0]+fOut[i][1]*fOut[i][1]));
pwrFourier->SetBinError(i+1, 0.0);
}
@@ -326,7 +344,13 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
snprintf(name, sizeof(name), "%s_Fourier_Phase", fData->GetName());
snprintf(title, sizeof(title), "%s_Fourier_Phase", fData->GetTitle());
TH1F* phaseFourier = new TH1F(name, title, fNoOfBins/2, -fResolution/2.0, (Double_t)fNoOfBins/2.0*fResolution+fResolution/2.0);
UInt_t noOfFourierBins = 0;
if (fNoOfBins % 2 == 0)
noOfFourierBins = fNoOfBins/2;
else
noOfFourierBins = (fNoOfBins+1)/2;
TH1F* phaseFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0);
if (phaseFourier == 0) {
fValid = false;
cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the phase part of the Fourier transform." << endl;
@@ -335,7 +359,7 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
// fill phaseFourier vector
Double_t value = 0.0;
for (UInt_t i=0; i<fNoOfBins/2; i++) {
for (UInt_t i=0; i<noOfFourierBins; i++) {
// calculate the phase
if (fOut[i][0] == 0) {
if (fOut[i][1] >= 0.0)
@@ -380,14 +404,23 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
}
}
// 2nd fill fIn
Int_t val = static_cast<Int_t>(fStartTime/fTimeResolution) + t0bin;
Int_t ival = static_cast<Int_t>(fStartTime/fTimeResolution) + t0bin;
UInt_t start = 0;
if (val >= 0) {
start = static_cast<UInt_t>(static_cast<Int_t>(fStartTime/fTimeResolution) + t0bin);
if (ival >= 0) {
start = static_cast<UInt_t>(ival);
}
Double_t mean = 0.0;
if (fDCCorrected) {
for (UInt_t i=0; i<fNoOfData; i++) {
mean += fData->GetBinContent(i+start);
}
mean /= (Double_t)fNoOfData;
}
// 2nd fill fIn
for (UInt_t i=0; i<fNoOfData; i++) {
fIn[i][0] = fData->GetBinContent(i+start);
fIn[i][0] = fData->GetBinContent(i+start) - mean;
fIn[i][1] = 0.0;
}
for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
src/classes/PFunction.cpp
+2 -1
View File
@@ -36,6 +36,7 @@ using namespace std;
#include <boost/algorithm/string/trim.hpp> // for stripping leading whitespace in std::string
#include "PMusr.h"
#include "PFunction.h"
//--------------------------------------------------------------------------
@@ -188,7 +189,7 @@ void PFunction::FillFuncEvalTree(iter_t const& i, PFuncTreeNode &node)
node.fDvalue = 3.14159265358979323846; // keep the value
} else if (i->value.id() == PFunctionGrammar::constGammaMuID) { // handle constant gamma_mu
node.fID = PFunctionGrammar::constGammaMuID; // keep the ID
node.fDvalue = 0.0135538817; // keep the value
node.fDvalue = GAMMA_BAR_MUON; // keep the value
} else if (i->value.id() == PFunctionGrammar::parameterID) { // handle parameter number
str = string(i->value.begin(), i->value.end()); // get string
boost::algorithm::trim(str);
src/classes/PMsrHandler.cpp
+24 -2
View File
@@ -935,6 +935,12 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout << endl;
} else if (sstr.BeginsWith("fourier_power")) {
fout << "fourier_power " << fFourier.fFourierPower << endl;
} else if (sstr.BeginsWith("dc-corrected")) {
fout << "dc-corrected ";
if (fFourier.fDCCorrected == true)
fout << "true" << endl;
else
fout << "false" << endl;
} else if (sstr.BeginsWith("apodization")) {
fout << "apodization ";
if (fFourier.fApodization == FOURIER_APOD_NONE) {
@@ -3156,6 +3162,7 @@ void PMsrHandler::InitFourierParameterStructure(PMsrFourierStructure &fourier)
fourier.fFourierBlockPresent = false; // fourier block present
fourier.fUnits = FOURIER_UNIT_NOT_GIVEN; // fourier untis, default: NOT GIVEN
fourier.fFourierPower = -1; // zero padding, default: -1 = NOT GIVEN
fourier.fDCCorrected = false; // dc-corrected FFT, default: false
fourier.fApodization = FOURIER_APOD_NOT_GIVEN; // apodization, default: NOT GIVEN
fourier.fPlotTag = FOURIER_PLOT_NOT_GIVEN; // initial plot tag, default: NOT GIVEN
fourier.fPhaseParamNo = 0; // initial parameter no = 0 means not a parameter
@@ -3187,8 +3194,6 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines)
if (lines.empty()) // no fourier block present
return true;
//cout << endl << ">> in PMsrHandler::HandleFourierEntry, Fourier block present ...";
PMsrFourierStructure fourier;
InitFourierParameterStructure(fourier);
@@ -3251,6 +3256,22 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines)
continue;
}
}
} else if (iter->fLine.BeginsWith("dc-corrected", TString::kIgnoreCase)) { // dc-corrected
if (tokens->GetEntries() < 2) { // dc-corrected tag is missing
error = true;
continue;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (!str.CompareTo("true", TString::kIgnoreCase) || !str.CompareTo("1")) {
fourier.fDCCorrected = true;
} else if (!str.CompareTo("false", TString::kIgnoreCase) || !str.CompareTo("0")) {
fourier.fDCCorrected = false;
} else { // unrecognized dc-corrected tag
error = true;
continue;
}
}
} else if (iter->fLine.BeginsWith("apodization", TString::kIgnoreCase)) { // apodization
if (tokens->GetEntries() < 2) { // apodization tag is missing
error = true;
@@ -3388,6 +3409,7 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines)
cerr << endl << "[fourier_power n # n is a number such that zero padding up to 2^n will be used]";
cerr << endl << " n=0 means no zero padding";
cerr << endl << " 0 <= n <= 20 are allowed values";
cerr << endl << "[dc-corrected true | false]";
cerr << endl << "[apodization none | weak | medium | strong]";
cerr << endl << "[plot real | imag | real_and_imag | power | phase]";
cerr << endl << "[phase value]";
src/classes/PMusrCanvas.cpp
+4 -3
View File
@@ -369,6 +369,7 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler)
if (fMsrHandler->GetMsrFourierList()->fFourierPower != -1) {
fFourier.fFourierPower = fMsrHandler->GetMsrFourierList()->fFourierPower;
}
fFourier.fDCCorrected = fMsrHandler->GetMsrFourierList()->fDCCorrected;
if (fMsrHandler->GetMsrFourierList()->fApodization != FOURIER_APOD_NOT_GIVEN) {
fFourier.fApodization = fMsrHandler->GetMsrFourierList()->fApodization;
}
@@ -2326,7 +2327,7 @@ void PMusrCanvas::HandleFourier()
double endTime = fHistoFrame->GetBinCenter(bin);
for (UInt_t i=0; i<fData.size(); i++) {
// calculate fourier transform of the data
PFourier fourierData(fData[i].data, fFourier.fUnits, startTime, endTime, fFourier.fFourierPower);
PFourier fourierData(fData[i].data, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData.IsValid()) {
cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier data ..." << endl;
return;
@@ -2366,7 +2367,7 @@ void PMusrCanvas::HandleFourier()
// calculate fourier transform of the theory
Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3;
PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, powerPad);
PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad);
if (!fourierTheory.IsValid()) {
cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier theory ..." << endl;
return;
@@ -2486,7 +2487,7 @@ void PMusrCanvas::HandleDifferenceFourier()
for (UInt_t i=0; i<fData.size(); i++) {
// calculate fourier transform of the data
PFourier fourierData(fData[i].diff, fFourier.fUnits, startTime, endTime, fFourier.fFourierPower);
PFourier fourierData(fData[i].diff, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower);
if (!fourierData.IsValid()) {
cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier diff ..." << endl;
return;
src/classes/PRunAsymmetry.cpp
+2 -2
View File
@@ -1604,10 +1604,10 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
gammaRRF = 1.0;
break;
case RRF_UNIT_G:
gammaRRF = 0.0135538817*TMath::TwoPi();
gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi();
break;
case RRF_UNIT_T:
gammaRRF = 0.0135538817*TMath::TwoPi()*1.0e4;
gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi()*1.0e4;
break;
default:
gammaRRF = TMath::TwoPi();
src/classes/PRunSingleHisto.cpp
+2 -3
View File
@@ -866,7 +866,6 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
// write these times back into the data structure. This way it is available when writting the log-file
fRunInfo->SetFitRange(fFitStartTime, 0);
fRunInfo->SetFitRange(fFitEndTime, 1);
cout << "debug> fit: " << fRunInfo->GetFitRange(0) << ", " << fRunInfo->GetFitRange(1) << end;
}
// check how the background shall be handled
@@ -1302,10 +1301,10 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
gammaRRF = 1.0;
break;
case RRF_UNIT_G:
gammaRRF = 0.0135538817*TMath::TwoPi();
gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi();
break;
case RRF_UNIT_T:
gammaRRF = 0.0135538817*TMath::TwoPi()*1.0e4;
gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi()*1.0e4;
break;
default:
gammaRRF = TMath::TwoPi();