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

back-integrate DKS Fourier after DKS upgrade and modularization. Not yet tested other than compilation.

Browse Source
This commit is contained in:
suter_a 2017-05-29 17:54:57 +02:00
parent c2e362d406
commit 2b7a391bbb
2 changed files with 3 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
src/classes/PFourier.cpp
View File

@ -8,7 +8,7 @@
***************************************************************************/ ***************************************************************************/
/*************************************************************************** /***************************************************************************
* Copyright (C) 2007-2016 by Andreas Suter * * Copyright (C) 2007-2017 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 *
@ -298,12 +298,10 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
fIn = 0; fIn = 0;
fOut = 0; fOut = 0;
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
fInDKS = 0; fInDKS = 0;
fOutDKS = 0; fOutDKS = 0;
#endif #endif
*/
SetUseFFTW(useFFTW); SetUseFFTW(useFFTW);
@ -365,13 +363,11 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins); fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
fOut = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins); fOut = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfBins);
} else { // try DKS } else { // try DKS
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
fInDKS = new double[fNoOfBins]; fInDKS = new double[fNoOfBins];
unsigned int size=fNoOfBins/2+1; unsigned int size=fNoOfBins/2+1;
fOutDKS = new complex<double>[size]; fOutDKS = new complex<double>[size];
#endif #endif
*/
} }
// check if memory allocation has been successful // check if memory allocation has been successful
@ -381,17 +377,12 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
return; return;
} }
} else { // try DKS } else { // try DKS
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
if ((fInDKS == 0) || (fOutDKS == 0)) { if ((fInDKS == 0) || (fOutDKS == 0)) {
fValid = false; fValid = false;
return; return;
} }
#else
fValid = false;
return;
#endif #endif
*/
fValid = false; fValid = false;
return; return;
} }
@ -405,7 +396,6 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
fValid = false; fValid = false;
} }
} else { // try DKS } else { // try DKS
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
// init DKSBase // init DKSBase
fDks.setAPI("Cuda", 4); fDks.setAPI("Cuda", 4);
@ -430,10 +420,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
} }
if ((fReal_ptr==0) || (fComp_ptr==0)) if ((fReal_ptr==0) || (fComp_ptr==0))
fValid = false; fValid = false;
#else
fValid = false;
#endif #endif
*/
fValid = false; fValid = false;
} }
} }
@ -454,7 +441,6 @@ PFourier::~PFourier()
if (fOut) if (fOut)
fftw_free(fOut); fftw_free(fOut);
} else { } else {
/*//as
#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);
@ -465,7 +451,6 @@ PFourier::~PFourier()
if (fOut) if (fOut)
delete [] fOutDKS; delete [] fOutDKS;
#endif #endif
*/
} }
} }
@ -488,7 +473,6 @@ void PFourier::Transform(UInt_t apodizationTag)
if (fUseFFTW) { if (fUseFFTW) {
fftw_execute(fFFTwPlan); fftw_execute(fFFTwPlan);
} else { } else {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
int dimsize[3] = {(int)fNoOfBins, 1, 1}; int dimsize[3] = {(int)fNoOfBins, 1, 1};
int status=0, size=fNoOfBins/2+1; int status=0, size=fNoOfBins/2+1;
@ -498,11 +482,7 @@ void PFourier::Transform(UInt_t apodizationTag)
status = fDks.callR2CFFT(fReal_ptr, fComp_ptr, 1, dimsize); status = fDks.callR2CFFT(fReal_ptr, fComp_ptr, 1, dimsize);
// read data from accelerator // read data from accelerator
status = fDks.readData< complex<double> >(fComp_ptr, fOutDKS, size); status = fDks.readData< complex<double> >(fComp_ptr, fOutDKS, size);
#else
fValid = false;
return;
#endif #endif
*/
fValid = false; fValid = false;
return; return;
} }
@ -527,7 +507,6 @@ void PFourier::Transform(UInt_t apodizationTag)
fOut[i][1] = im; fOut[i][1] = im;
} }
} else { // try DKS } else { // try DKS
/*//as
UInt_t size=fNoOfBins/2+1; UInt_t size=fNoOfBins/2+1;
for (UInt_t i=0; i<size; i++) { for (UInt_t i=0; i<size; i++) {
phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime; phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime;
@ -537,7 +516,6 @@ void PFourier::Transform(UInt_t apodizationTag)
fOutDKS[i] = complex<double>(re, im); fOutDKS[i] = complex<double>(re, im);
#endif #endif
} }
*/
} }
} }
@ -609,15 +587,11 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
} }
} else { } else {
for (UInt_t i=0; i<noOfFourierBins; i++) { for (UInt_t i=0; i<noOfFourierBins; i++) {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
realFourier->SetBinContent(i+1, scale*fOutDKS[i].real()); realFourier->SetBinContent(i+1, scale*fOutDKS[i].real());
#else #else
realFourier->SetBinContent(i+1, PMUSR_UNDEFINED); realFourier->SetBinContent(i+1, PMUSR_UNDEFINED);
#endif #endif
*/
realFourier->SetBinContent(i+1, PMUSR_UNDEFINED); //as needs to be removed
realFourier->SetBinError(i+1, 0.0);
} }
} }
return realFourier; return realFourier;
@ -761,15 +735,11 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
} }
} else { } else {
for (UInt_t i=0; i<noOfFourierBins; i++) { for (UInt_t i=0; i<noOfFourierBins; i++) {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
imaginaryFourier->SetBinContent(i+1, scale*fOutDKS[i].imag()); imaginaryFourier->SetBinContent(i+1, scale*fOutDKS[i].imag());
#else #else
imaginaryFourier->SetBinContent(i+1, PMUSR_UNDEFINED); imaginaryFourier->SetBinContent(i+1, PMUSR_UNDEFINED);
#endif #endif
*/
imaginaryFourier->SetBinContent(i+1, PMUSR_UNDEFINED); //as needs to be removed
imaginaryFourier->SetBinError(i+1, 0.0);
} }
} }
return imaginaryFourier; return imaginaryFourier;
@ -812,15 +782,11 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
} }
} else { } else {
for (UInt_t i=0; i<noOfFourierBins; i++) { for (UInt_t i=0; i<noOfFourierBins; i++) {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
pwrFourier->SetBinContent(i+1, scale*sqrt(fOutDKS[i].real()*fOutDKS[i].real()+fOutDKS[i].imag()*fOutDKS[i].imag())); pwrFourier->SetBinContent(i+1, scale*sqrt(fOutDKS[i].real()*fOutDKS[i].real()+fOutDKS[i].imag()*fOutDKS[i].imag()));
#else #else
pwrFourier->SetBinContent(i+1, PMUSR_UNDEFINED); pwrFourier->SetBinContent(i+1, PMUSR_UNDEFINED);
#endif #endif
*/
pwrFourier->SetBinContent(i+1, PMUSR_UNDEFINED); //as needs to be removed
pwrFourier->SetBinError(i+1, 0.0);
} }
} }
return pwrFourier; return pwrFourier;
@ -863,7 +829,6 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
re = fOut[i][0]; re = fOut[i][0];
im = fOut[i][1]; im = fOut[i][1];
} else { } else {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
re = fOutDKS[i].real(); re = fOutDKS[i].real();
im = fOutDKS[i].imag(); im = fOutDKS[i].imag();
@ -871,9 +836,6 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
re = 1.0; re = 1.0;
im = 0.0; im = 0.0;
#endif #endif
*/
re = 1.0; //as needs to be removed
im = 0.0; //as needs to be removed
} }
// calculate the phase // calculate the phase
if (re == 0) { if (re == 0) {
@ -946,18 +908,14 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
} }
} else { } else {
for (UInt_t i=0; i<fNoOfData-t0bin; i++) { for (UInt_t i=0; i<fNoOfData-t0bin; i++) {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
fInDKS[i] = fData->GetBinContent(i+t0bin) - mean; fInDKS[i] = fData->GetBinContent(i+t0bin) - mean;
#endif #endif
*/
} }
for (UInt_t i=fNoOfData-t0bin; i<fNoOfBins; i++) { for (UInt_t i=fNoOfData-t0bin; i<fNoOfBins; i++) {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
fInDKS[i] = 0.0; fInDKS[i] = 0.0;
#endif #endif
*/
} }
} }
@ -1021,11 +979,9 @@ void PFourier::ApodizeData(Int_t apodizationTag) {
if (fUseFFTW) { if (fUseFFTW) {
fIn[i][0] *= q; fIn[i][0] *= q;
} else { } else {
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
fInDKS[i] *= q; fInDKS[i] *= q;
#endif #endif
*/
} }
} }
} }

11
src/include/PFourier.h
View File

@ -34,11 +34,6 @@
#include "config.h" #include "config.h"
#endif #endif
//as the next 3 lines will need to be deleted as soon as DKS Fourier is back in place
#include <vector>
using namespace std;
#include "fftw3.h"
/*//as
#ifndef HAVE_DKS #ifndef HAVE_DKS
#include <vector> #include <vector>
using namespace std; using namespace std;
@ -47,8 +42,8 @@ using namespace std;
#include <complex> #include <complex>
using namespace std; using namespace std;
#include "DKSBase.h" #include "DKSBase.h"
#include "DKSFFT.h"
#endif #endif
*/
#include "Minuit2/FCNBase.h" #include "Minuit2/FCNBase.h"
@ -159,15 +154,13 @@ class PFourier
fftw_complex *fIn; ///< real part of the Fourier transform fftw_complex *fIn; ///< real part of the Fourier transform
fftw_complex *fOut; ///< imaginary part of the Fourier transform fftw_complex *fOut; ///< imaginary part of the Fourier transform
/*//as
#ifdef HAVE_DKS #ifdef HAVE_DKS
double *fInDKS; ///< real part of the Fourier transform double *fInDKS; ///< real part of the Fourier transform
complex<double> *fOutDKS; ///< imaginary part of the Fourier transform complex<double> *fOutDKS; ///< imaginary part of the Fourier transform
DKSBase fDks; ///< Dynamic Kernel Scheduler DKSFFT fDks; ///< Dynamic Kernel Scheduler
void *fReal_ptr; ///< real part of the Fourier on accelartor void *fReal_ptr; ///< real part of the Fourier on accelartor
void *fComp_ptr; ///< imaginary part of the Fourier on the acclerator void *fComp_ptr; ///< imaginary part of the Fourier on the acclerator
#endif #endif
*/
virtual void PrepareFFTwInputData(UInt_t apodizationTag); virtual void PrepareFFTwInputData(UInt_t apodizationTag);
virtual void ApodizeData(Int_t apodizationTag); virtual void ApodizeData(Int_t apodizationTag);