From 18a3055e7c66359b7e8ce4c3f416d51eb6b4af62 Mon Sep 17 00:00:00 2001 From: nemu Date: Tue, 21 Apr 2009 06:43:25 +0000 Subject: [PATCH] fixed a filling bug in PFourier --- src/classes/PFourier.cpp | 52 +++++++++++++++++++++++++++++++------ src/classes/PMusrCanvas.cpp | 19 +++++++++++--- 2 files changed, 59 insertions(+), 12 deletions(-) diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index aa6e189a..643b7179 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -63,6 +63,14 @@ PFourier::PFourier(TH1F *data, int unitTag, double startTime, double endTime, un return; } +//cout << endl << "PFourier::PFourier: fData = " << fData; +/* +for (unsigned int i=0; i<10; i++) { + cout << endl << "PFourier::PFourier: i=" << i << ", t=" << fData->GetBinCenter(i) << ", value=" << fData->GetBinContent(i); +} +*/ +//cout << endl; + fValid = true; fIn = 0; fOut = 0; @@ -71,7 +79,7 @@ PFourier::PFourier(TH1F *data, int unitTag, double startTime, double endTime, un // calculate time resolution in (us) fTimeResolution = fData->GetBinWidth(1); -cout << endl << ">> fTimeResolution = " << fTimeResolution; +//cout << endl << ">> fTimeResolution = " << fTimeResolution; // if endTime == 0 set it to the last time slot if (fEndTime == 0.0) { @@ -105,7 +113,7 @@ cout << endl << ">> fTimeResolution = " << fTimeResolution; fNoOfBins = fNoOfData; } -cout << endl << ">> fNoOfBins = " << fNoOfBins; +//cout << endl << ">> fNoOfBins = " << fNoOfBins; // calculate fourier resolution double resolution = 1.0/(fTimeResolution*fNoOfBins); // in MHz @@ -125,7 +133,7 @@ cout << endl << ">> fNoOfBins = " << fNoOfBins; break; } -cout << endl << ">> fResolution = " << fResolution; +//cout << endl << ">> fResolution = " << fResolution; // allocate necessary memory fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins); @@ -140,6 +148,7 @@ cout << endl << ">> fResolution = " << fResolution; return; } +//cout << endl << "PFourier::PFourier: fNoOfBins=" << fNoOfBins << endl; fFFTwPlan = fftw_plan_dft_1d(fNoOfBins, fIn, fOut, FFTW_FORWARD, FFTW_ESTIMATE); if (!fFFTwPlan) { @@ -159,7 +168,7 @@ cout << endl << ">> fResolution = " << fResolution; */ PFourier::~PFourier() { -cout << endl << "in ~PFourier() ..." << endl; +//cout << endl << "in ~PFourier() ..." << endl; if (fValid) { fftw_destroy_plan(fFFTwPlan); fftw_free(fIn); @@ -180,11 +189,27 @@ void PFourier::Transform(unsigned int apodizationTag) if (!fValid) return; -cout << endl << ">> PFourier::Transform ..." << endl; +//cout << endl << ">> PFourier::Transform (apodizationTag=" << apodizationTag << ") ..." << endl; PrepareFFTwInputData(apodizationTag); +/* +if (fNoOfBins < 200) { +for (unsigned int i=0; i> PFourier::PrepareFFTwInputData: " << i << ": fIn[i][0] = " << fIn[i][0]; +} +cout << endl; +} +*/ + fftw_execute(fFFTwPlan); + +/* +for (unsigned int i=fNoOfBins-10; i> PFourier::PrepareFFTwInputData: " << i << ": fOut[i][0] = " << fOut[i][0]; +} +cout << endl; +*/ } //-------------------------------------------------------------------------- @@ -219,7 +244,7 @@ TH1F* PFourier::GetRealFourier(const double scale) realFourier->SetBinContent(i+1, scale*fOut[i][0]); realFourier->SetBinError(i+1, 0.0); } -cout << endl << ">> realFourier = " << realFourier; +//cout << endl << ">> realFourier = " << realFourier; return realFourier; } @@ -359,7 +384,9 @@ void PFourier::PrepareFFTwInputData(unsigned int apodizationTag) { // 1st find t==0. fData start at times t<0!! int t0bin = -1; - for (int i=0; iGetNbinsX(); i++) { +//cout << ">> PFourier::PrepareFFTwInputData: fData=" << fData << ", fData->GetNbinsX() = " << fData->GetNbinsX(); + for (int i=0; iGetNbinsX(); i++) { +//if (i<20) cout << endl << ">> PFourier::PrepareFFTwInputData: i=" << i << ", fData->GetBinCenter(i)=" << fData->GetBinCenter(i); if (fData->GetBinCenter(i) >= 0.0) { t0bin = i; break; @@ -372,13 +399,22 @@ void PFourier::PrepareFFTwInputData(unsigned int apodizationTag) fIn[i][0] = fData->GetBinContent(i+start+1); fIn[i][1] = 0.0; } - for (unsigned int i=fNoOfData; i> PFourier::PrepareFFTwInputData: t0bin = " << t0bin << ", start = " << start << endl; // 3rd apodize data (if wished) ApodizeData(apodizationTag); + +/* +cout << ">> PFourier::PrepareFFTwInputData: fNoOfData = " << fNoOfData << ", fNoOfBins = " << fNoOfBins << endl; +for (unsigned int i=0; i<10; i++) { + cout << endl << ">> PFourier::PrepareFFTwInputData: " << i << ": fIn[i][0] = " << fIn[i][0]; +} +cout << endl; +*/ } //-------------------------------------------------------------------------- diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index b60a51a9..a0c18d38 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1568,12 +1568,16 @@ void PMusrCanvas::HandleFourier() // check if fourier needs to be calculated if (fData[0].dataFourierRe == 0) { +//cout << endl << ">> Recalculate Fourier ----------------------------------------"; +//cout << endl << ">> fData[0].data = " << fData[0].data; int bin; bin = fData[0].data->GetXaxis()->GetFirst(); +//cout << endl << ">> start bin = " << bin; double startTime = fData[0].data->GetBinCenter(bin); +//cout << endl << ">> start time = " << startTime; bin = fData[0].data->GetXaxis()->GetLast(); double endTime = fData[0].data->GetBinCenter(bin); -//cout << endl << ">> startTime = " << startTime << ", endTime = " << endTime << endl; +//cout << endl << ">> Fourier: startTime = " << startTime << ", endTime = " << endTime; for (unsigned int i=0; i> Fourier: i=" << i; +for (unsigned j=0; j<10; j++) { + cout << endl << ">> Fourier: " << j << ", data = " << fData[i].data->GetBinContent(j) << ", fourier.power = " << fData[i].dataFourierPwr->GetBinContent(j); +} +cout << endl; +*/ // set marker and line color fData[i].dataFourierRe->SetMarkerColor(fData[i].data->GetMarkerColor()); @@ -1618,7 +1629,7 @@ void PMusrCanvas::HandleFourier() // calculate fourier transform of the theory int powerPad = (int)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3; -cout << endl << ">> powerPad = " << powerPad; +//cout << endl << ">> powerPad = " << powerPad; PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, powerPad); if (!fourierTheory.IsValid()) { cout << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier theory ..." << endl; @@ -1626,7 +1637,7 @@ cout << endl << ">> powerPad = " << powerPad; } fourierTheory.Transform(fFourier.fApodization); scale = sqrt(fData[0].theory->GetBinWidth(1)/(endTime-startTime)*fData[0].theory->GetBinWidth(1)/fData[0].data->GetBinWidth(1)); -cout << endl << ">> theory scale = " << scale << ", data.res/theory.res = " << fData[0].theory->GetBinWidth(1)/fData[0].data->GetBinWidth(1); +//cout << endl << ">> theory scale = " << scale << ", data.res/theory.res = " << fData[0].theory->GetBinWidth(1)/fData[0].data->GetBinWidth(1); // get real part of the data fData[i].theoryFourierRe = fourierTheory.GetRealFourier(scale); //cout << endl << ">> i: " << i << ", fData[i].dataFourierRe = " << fData[i].dataFourierRe; @@ -2313,7 +2324,7 @@ void PMusrCanvas::PlotFourier() xAxisTitle = TString("??"); } - // plot data + // plot fourier data double min, max, binContent; switch (fCurrentPlotView) { case PV_FOURIER_REAL: