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

added apodization

Browse Source
This commit is contained in:
nemu
2008-09-04 07:13:04 +00:00
parent 0277420671
commit e51acbe0fb
3 changed files with 148 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

116
src/tests/fourier/PMusrFourier.cpp
View File

@ -75,6 +75,9 @@ PMusrFourier::PMusrFourier(int dataType, PDoubleVector &data, double timeResolut
fIn = 0;
fOut = 0;
fApodization = F_APODIZATION_NONE;
fFilter = F_FILTER_NONE;
// if endTime == 0 set is to the last time slot
if (fEndTime == 0.0)
fEndTime = (fData.size()-1)*fTimeResolution;
@ -106,7 +109,7 @@ PMusrFourier::PMusrFourier(int dataType, PDoubleVector &data, double timeResolut
}
}
cout << endl << "dB = " << 1.0/(F_GAMMA_BAR_MUON * (fEndTime-fStartTime)) << " (G), Bmax = " << 1.0/(F_GAMMA_BAR_MUON*fTimeResolution) << " (G)" << endl;
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
@ -175,17 +178,18 @@ cout << endl << "in ~PMusrFourier() ..." << endl;
/**
* <p>
*
* \param apodization
* \param filter
* \param estimateN0AndBkg
* \param apodizationTag 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
* \param filterTag 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
*/
void PMusrFourier::Transform(int apodization, int filter)
void PMusrFourier::Transform(unsigned int apodizationTag, unsigned int filterTag)
{
if (!fValid)
return;
if (fDataType == F_SINGLE_HISTO) {
PrepareSingleHistoFFTwInputData();
PrepareSingleHistoFFTwInputData(apodizationTag, filterTag);
} else {
PrepareFFTwInputData(apodizationTag, filterTag);
}
// for test only
@ -212,7 +216,7 @@ for (int p=-180; p<180; p++) {
}
fftw_execute(fFFTwPlan);
if (p==8) {
if (p == 7) {
for (unsigned int j=0; j<fNoOfBins/2; j++) {
re.SetBinContent(j+1, fOut[j][0]);
im.SetBinContent(j+1, fOut[j][1]);
@ -224,6 +228,9 @@ for (int p=-180; p<180; p++) {
for (unsigned int i=0; i<fNoOfBins/2; i++) {
sum += fOut[i][1];
}
if (p == 7) {
cout << endl << ">> sum = " << sum << endl;
}
sumHist.SetBinContent(p+181, fabs(sum));
}
TFile f("test_out.root", "RECREATE");
@ -406,7 +413,7 @@ cout << endl << ">> N0/per bin=" << A/(PMUON_LIFETIME*1000.0)*fTimeResolution <<
* <p>
*
*/
void PMusrFourier::PrepareSingleHistoFFTwInputData()
void PMusrFourier::PrepareSingleHistoFFTwInputData(unsigned int apodizationTag, unsigned int filterTag)
{
// 1st fill fIn
unsigned int start = (unsigned int)(fStartTime/fTimeResolution);
@ -431,6 +438,9 @@ void PMusrFourier::PrepareSingleHistoFFTwInputData()
for (unsigned int i=0; i<fNoOfData; i++)
fIn[i][0] -= fN0;
// 5th apodize data (if wished)
ApodizeData(apodizationTag);
// has to be removed after testing
TH1F test_in("test_in", "test_in", fNoOfBins,
fStartTime - fTimeResolution/2.0, fNoOfBins*fTimeResolution + fTimeResolution/2.0);
@ -440,3 +450,93 @@ TFile f("test_in.root", "RECREATE");
test_in.Write();
f.Close();
}
//--------------------------------------------------------------------------
// PrepareFFTwInputData
//--------------------------------------------------------------------------
/**
* <p>
*
*/
void PMusrFourier::PrepareFFTwInputData(unsigned int apodizationTag, unsigned int filterTag)
{
// 1st fill fIn
unsigned int start = (unsigned int)(fStartTime/fTimeResolution);
for (unsigned int i=0; i<fNoOfData; i++) {
fIn[i][0] = fDataRebinned[i+start];
fIn[i][1] = 0.0;
}
for (unsigned int i=fNoOfData; i<fNoOfBins; i++) {
fIn[i][0] = 0.0;
fIn[i][1] = 0.0;
}
// 2nd apodize data (if wished)
ApodizeData(apodizationTag);
// has to be removed after testing
TH1F test_in("test_in", "test_in", fNoOfBins,
fStartTime - fTimeResolution/2.0, fNoOfBins*fTimeResolution + fTimeResolution/2.0);
for (unsigned int i=0; i<fNoOfBins; i++)
test_in.SetBinContent(i, fIn[i][0]);
TFile f("test_in.root", "RECREATE");
test_in.Write();
f.Close();
}
//--------------------------------------------------------------------------
// ApodizeData
//--------------------------------------------------------------------------
/**
* <p>
*
* \param apodizationTag
*/
void PMusrFourier::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:
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;
case F_APODIZATION_USER:
cout << endl << ">> User Apodization not yet implemented, sorry ..." << endl;
break;
default:
cout << endl << ">> **ERROR** User Apodization tag " << apodizationTag << " unknown, sorry ..." << endl;
break;
}
double q;
for (unsigned int i=0; i<fNoOfData; i++) {
q = c[0];
for (unsigned int j=1; j<5; j++) {
q += c[j] * pow((double)i/(double)fNoOfData, 2.0*(double)j);
}
fIn[i][0] *= q;
}
}

17
src/tests/fourier/PMusrFourier.h
View File

@ -39,8 +39,9 @@ using namespace std;
#define F_ESTIMATE_N0_AND_BKG true
#define F_SINGLE_HISTO 0
#define F_ASYMMETRY 1
#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
@ -87,8 +88,9 @@ class PMusrFourier
virtual void SetN0(double n0) { fN0 = n0; }
virtual void SetBkg(double bkg) { fBkg = bkg; }
virtual void SetPhase(double phase) { fPhase = phase; }
virtual void Transform(int apodization = -1, int filter = -1);
virtual void Transform(unsigned int apodizationTag = 0, unsigned int filter = 0);
virtual double GetFieldResolution() { return fFieldResolution; }
virtual double GetPhaseCorrection() { return fPhaseCorrection; }
@ -103,11 +105,12 @@ class PMusrFourier
bool fValid;
int fDataType; ///< 0=Single Histo, 1=Asymmetry
int fUseApodization; ///< 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
int fUseFilter; ///< 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
int fApodization; ///< 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
int fFilter; ///< 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
double fN0;
double fBkg;
double fPhase;
double fTimeResolution;
double fStartTime;
@ -126,7 +129,9 @@ class PMusrFourier
fftw_complex *fIn;
fftw_complex *fOut;
virtual void PrepareSingleHistoFFTwInputData();
virtual void PrepareSingleHistoFFTwInputData(unsigned int apodizationTag, unsigned int filterTag);
virtual void PrepareFFTwInputData(unsigned int apodizationTag, unsigned int filterTag);
virtual void ApodizeData(int apodizationTag);
virtual void EstimateN0AndBkg();
};

31
src/tests/fourier/fourier.cpp
View File

@ -159,11 +159,38 @@ cout << endl << "#bins=" << histo->GetNbinsX();
cin >> zeroPaddingPower;
}
bool estimateN0AndBkg = true;
double N0, bkg, phase;
cout << endl << ">> Do you want to give N0, Bkg, and phase explicitly (y/n)? ";
cin >> answer;
if (strstr(answer, "y")) {
estimateN0AndBkg = false;
cout << endl << ">> N0 = ";
cin >> N0;
cout << endl << ">> bkg = ";
cin >> bkg;
cout << endl << ">> phase = ";
cin >> phase;
}
PMusrFourier fourier(F_SINGLE_HISTO, data, timeResolution, startTime, endTime,
rebin, zeroPaddingPower, F_ESTIMATE_N0_AND_BKG);
rebin, zeroPaddingPower, estimateN0AndBkg);
if (fourier.IsValid()) {
fourier.Transform(F_APODIZATION_NONE, F_FILTER_NONE);
if (!estimateN0AndBkg) {
fourier.SetN0(N0);
fourier.SetBkg(bkg);
fourier.SetPhase(phase);
}
cout << endl << ">> Do you wish to apodize your data (y/n)? ";
cin >> answer;
unsigned int apodizationTag=0;
if (strstr(answer, "y")) {
cout << endl << ">> apodization (1=weak, 2=medium, 3=strong, 4=user) = ";
cin >> apodizationTag;
}
fourier.Transform(apodizationTag, F_FILTER_NONE);
} else {
cout << endl << "**ERROR** something went wrong when invoking the PMusrFourier object ...";
cout << endl << endl;