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

fixed wrong Fourier vector size for DKS.

Browse Source
This commit is contained in:
suter_a 2015-04-08 17:24:24 +02:00
parent 75578f1977
commit a393cb9ec5
2 changed files with 8 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
src/classes/PFourier.cpp
View File

@ -8,7 +8,7 @@
***************************************************************************/ ***************************************************************************/
/*************************************************************************** /***************************************************************************
* Copyright (C) 2007-2014 by Andreas Suter * * Copyright (C) 2007-2015 by Andreas Suter *
* andreas.suter@psi.ch * * andreas.suter@psi.ch *
* * * *
* This program is free software; you can redistribute it and/or modify * * This program is free software; you can redistribute it and/or modify *
@ -99,6 +99,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
// calculate start and end bin // calculate start and end bin
fNoOfData = (UInt_t)((fEndTime-fStartTime)/fTimeResolution); fNoOfData = (UInt_t)((fEndTime-fStartTime)/fTimeResolution);
cout << "debug> fStartTime=" << fStartTime << ", fEndTime=" << fEndTime << ", fTimeResolution=" << fTimeResolution << ", fNoOfData=" << fNoOfData << endl;
// check if zero padding is whished // check if zero padding is whished
if (fZeroPaddingPower > 0) { if (fZeroPaddingPower > 0) {
@ -139,9 +140,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
} else { // try DKS } else { // try DKS
#ifdef HAVE_DKS #ifdef HAVE_DKS
fInDKS = new double[fNoOfBins]; fInDKS = new double[fNoOfBins];
unsigned int size=fNoOfBins/2; unsigned int size=fNoOfBins/2+1;
if (fNoOfBins % 2 == 1)
size++;
fOutDKS = new complex<double>[size]; fOutDKS = new complex<double>[size];
#endif #endif
} }
@ -182,9 +181,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
fDks.setupFFT(1, dimsize); fDks.setupFFT(1, dimsize);
// allocate memory on accelerator // allocate memory on accelerator
int ierr; int ierr;
unsigned int size=fNoOfBins/2; unsigned int size=fNoOfBins/2+1;
if (fNoOfBins % 2 == 1)
size++;
fReal_ptr = 0; fReal_ptr = 0;
fComp_ptr = 0; fComp_ptr = 0;
fReal_ptr = fDks.allocateMemory<double>(fNoOfBins, ierr); fReal_ptr = fDks.allocateMemory<double>(fNoOfBins, ierr);
@ -216,9 +213,7 @@ PFourier::~PFourier()
#ifdef HAVE_DKS #ifdef HAVE_DKS
// free accelerator memory // free accelerator memory
fDks.freeMemory<double>(fReal_ptr, (int)fNoOfBins); fDks.freeMemory<double>(fReal_ptr, (int)fNoOfBins);
int size = fNoOfBins/2; int size = fNoOfBins/2+1;
if (fNoOfBins % 2 == 1)
size++;
fDks.freeMemory< complex<double> >(fComp_ptr, size); fDks.freeMemory< complex<double> >(fComp_ptr, size);
if (fIn) if (fIn)
delete [] fInDKS; delete [] fInDKS;
@ -249,9 +244,7 @@ void PFourier::Transform(UInt_t apodizationTag)
} else { } else {
#ifdef HAVE_DKS #ifdef HAVE_DKS
int dimsize[3] = {fNoOfBins, 1, 1}; int dimsize[3] = {fNoOfBins, 1, 1};
int status=0, size=fNoOfBins/2; int status=0, size=fNoOfBins/2+1;
if (fNoOfBins % 2 == 1)
size++;
// write data to the accelerator // write data to the accelerator
status=fDks.writeData<double>(fReal_ptr, fInDKS, fNoOfBins); status=fDks.writeData<double>(fReal_ptr, fInDKS, fNoOfBins);
cout << "debug> status=" << status << ": write" << endl; cout << "debug> status=" << status << ": write" << endl;

3
src/tests/DKS_FourierTest/dks_fourierTest.cpp
View File

@ -78,7 +78,7 @@ int main(int argc, char *argv[])
snprintf(str, sizeof(str), "h%d", i); snprintf(str, sizeof(str), "h%d", i);
h = new TH1F(str, str, L+1, -dt/2, 10.0+dt/2); h = new TH1F(str, str, L+1, -dt/2, 10.0+dt/2);
for (unsigned int j=0; j<(unsigned int)h->GetNbinsX(); j++) { for (unsigned int j=0; j<(unsigned int)h->GetNbinsX(); j++) {
dval = exp(-0.5*dt*j)*cos(0.013554*5000.0*dt*j+6.2832/N*i); dval = exp(-0.25*pow(dt*j,2.0))*cos(0.013554*5000.0*dt*j+6.2832/N*i);
h->SetBinContent(j+1, dval); h->SetBinContent(j+1, dval);
} }
data.push_back(h); data.push_back(h);
@ -125,6 +125,7 @@ int main(int argc, char *argv[])
} }
TFile fout("dks_fourierTest.root", "recreate"); TFile fout("dks_fourierTest.root", "recreate");
for (unsigned int i=0; i<fourierPowerHistos.size(); i++) { for (unsigned int i=0; i<fourierPowerHistos.size(); i++) {
data[i]->Write();
fourierPowerHistos[i]->Write(); fourierPowerHistos[i]->Write();
} }
fout.Close(); fout.Close();