corrected size for noOfFourierBins, and made switching between FFTW and DKS more coherent

src/classes/PFourier.cpp

@@ -99,6 +99,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
 
   // calculate start and end bin
   fNoOfData = (UInt_t)((fEndTime-fStartTime)/fTimeResolution);
+  cout << "debug> fNoOfData=" << fNoOfData << endl;
 
   // check if zero padding is whished
   if (fZeroPaddingPower > 0) {
@@ -184,7 +185,15 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi
     fReal_ptr = 0;
     fComp_ptr = 0;
     fReal_ptr = fDks.allocateMemory<double>(fNoOfBins, ierr);
+    if (ierr > 0) {
+      cerr << ">> PFourier: **ERROR** Couldn't allocate memory for fReal_ptr." << endl;
+      fValid = false;
+    }
     fComp_ptr = fDks.allocateMemory< complex<double> >(size, ierr);
+    if (ierr > 0) {
+      cerr << ">> PFourier: **ERROR** Couldn't allocate memory for fComp_ptr." << endl;
+      fValid = false;
+    }
     if ((fReal_ptr==0) || (fComp_ptr==0))
       fValid = false;
 #else
@@ -245,7 +254,7 @@ void PFourier::Transform(UInt_t apodizationTag)
     int dimsize[3] = {fNoOfBins, 1, 1};
     int status=0, size=fNoOfBins/2+1;
     // write data to the accelerator
-    status=fDks.writeData<double>(fReal_ptr, fInDKS, fNoOfBins);
+    status = fDks.writeData<double>(fReal_ptr, fInDKS, fNoOfBins);
     // execute the FFT
     status = fDks.callR2CFFT(fReal_ptr, fComp_ptr, 1, dimsize);
     // read data from accelerator
@@ -276,7 +285,8 @@ void PFourier::Transform(UInt_t apodizationTag)
       fOut[i][1] = im;
     }
   } else { // try DKS
-    for (UInt_t i=0; i<fNoOfBins; i++) {
+    UInt_t size=fNoOfBins/2+1;
+    for (UInt_t i=0; i<size; i++) {
       phase = 2.0*PI/(fTimeResolution*fNoOfBins) * i * shiftTime;
 #ifdef HAVE_DKS
       re =  fOutDKS[i].real() * cos(phase) + fOutDKS[i].imag() * sin(phase);
@@ -314,11 +324,7 @@ void PFourier::SetUseFFTW(const Bool_t flag)
  */
 Double_t PFourier::GetMaxFreq()
 {
-  UInt_t noOfFourierBins = 0;
+  UInt_t noOfFourierBins = fNoOfBins/2+1;
-  if (fNoOfBins % 2 == 0)
-    noOfFourierBins = fNoOfBins/2;
-  else
-    noOfFourierBins = (fNoOfBins+1)/2;
 
   return fResolution*noOfFourierBins;
 }
@@ -343,11 +349,7 @@ 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());
 
-  UInt_t noOfFourierBins = 0;
+  UInt_t noOfFourierBins = fNoOfBins/2+1;
-  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) {
@@ -395,11 +397,7 @@ 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());
 
-  UInt_t noOfFourierBins = 0;
+  UInt_t noOfFourierBins = fNoOfBins/2+1;
-  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) {
@@ -447,11 +445,7 @@ 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());
 
-  UInt_t noOfFourierBins = 0;
+  UInt_t noOfFourierBins = fNoOfBins/2+1;
-  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) {
@@ -499,11 +493,7 @@ 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());
 
-  UInt_t noOfFourierBins = 0;
+  UInt_t noOfFourierBins = fNoOfBins/2+1;
-  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) {
@@ -564,45 +554,41 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
 void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
 {
   // 1st find t==0. fData start at times t<0!!
-  Int_t t0bin = -1;
+  UInt_t t0bin = -1;
   for (Int_t i=1; i<fData->GetNbinsX(); i++) {
     if (fData->GetBinCenter(i) >= 0.0) {
       t0bin = i;
       break;
     }
   }
-
+  cout << "debug> t0bin=" << t0bin << ", fNoOfData=" << fNoOfData << ", fNoOfBins=" << fNoOfBins << endl;
-  Int_t ival = static_cast<Int_t>(fStartTime/fTimeResolution) + t0bin;
-  UInt_t start = 0;
-  if (ival >= 0) {
-    start = static_cast<UInt_t>(ival);
-  }
 
   Double_t mean = 0.0;
   if (fDCCorrected) {
-    for (UInt_t i=0; i<fNoOfData; i++) {
+    for (UInt_t i=t0bin; i<fNoOfData; i++) {
-      mean += fData->GetBinContent(i+start);
+      mean += fData->GetBinContent(i);
     }
-    mean /= (Double_t)fNoOfData;
+    mean /= (Double_t)(fNoOfData-t0bin);
-  }
+    cout << "debug> mean=" << mean << endl;
+  }  
 
   // 2nd fill fIn
   if (fUseFFTW) {
-    for (UInt_t i=0; i<fNoOfData; i++) {
+    for (UInt_t i=0; i<fNoOfData-t0bin; i++) {
-      fIn[i][0] = fData->GetBinContent(i+start) - mean;
+      fIn[i][0] = fData->GetBinContent(i+t0bin) - mean;
       fIn[i][1] = 0.0;
     }
-    for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
+    for (UInt_t i=fNoOfData-t0bin; i<fNoOfBins; i++) {
       fIn[i][0] = 0.0;
       fIn[i][1] = 0.0;
     }
   } else {
-    for (UInt_t i=0; i<fNoOfData; i++) {
+    for (UInt_t i=0; i<fNoOfData-t0bin; i++) {
 #ifdef HAVE_DKS
-      fInDKS[i] = fData->GetBinContent(i+start) - mean;
+      fInDKS[i] = fData->GetBinContent(i+t0bin) - mean;
 #endif
     }
-    for (UInt_t i=fNoOfData; i<fNoOfBins; i++) {
+    for (UInt_t i=fNoOfData-t0bin; i<fNoOfBins; i++) {
 #ifdef HAVE_DKS
       fInDKS[i] = 0.0;
 #endif
@@ -625,6 +611,9 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
  */
 void PFourier::ApodizeData(Int_t apodizationTag) {
 
+  if (apodizationTag == F_APODIZATION_NONE)
+    return;
+
   const Double_t cweak[3]   = { 0.384093, -0.087577, 0.703484 };
   const Double_t cmedium[3] = { 0.152442, -0.136176, 0.983734 };
   const Double_t cstrong[3] = { 0.045335,  0.554883, 0.399782 };
@@ -635,9 +624,6 @@ void PFourier::ApodizeData(Int_t apodizationTag) {
   }
 
   switch (apodizationTag) {
-    case F_APODIZATION_NONE:
-      return;
-      break;
     case F_APODIZATION_WEAK:
       c[0] = cweak[0]+cweak[1]+cweak[2];
       c[1] = -(cweak[1]+2.0*cweak[2]);

src/musrFT.cpp

@@ -54,6 +54,10 @@ using namespace std;
 #include "PFourier.h"
 #include "PFourierCanvas.h"
 
+#define MFT_UNDEF -1
+#define MFT_DKS    0
+#define MFT_FFTW   1
+
 //----------------------------------------------------------------------------
 /**
  * <p>Structure keeping the command line options.
@@ -81,7 +85,7 @@ typedef struct {
   Int_t packing;                  ///< packing for rebinning the time histograms before Fourier transform.
   TString title;                  ///< title to be shown for the Fourier plot.
   Double_t lifetimecorrection;    ///< is == 0.0 for NO life time correction, otherwise it holds the fudge factor
-  Bool_t useDKS;                  ///< used DKS flag, true -> used DKS, false -> use FFTW
+  Int_t useDKS;                   ///< use DKS tag, -1 = undefined, 0 = used DKS, 1 = use FFTW
   Int_t timeout;                  ///< timeout in (sec) after which musrFT will terminate. if <= 0, no automatic termination will take place.
 } musrFT_startup_param;
 
@@ -182,7 +186,7 @@ void musrFT_init(musrFT_startup_param &startupParam)
   startupParam.packing = 1;
   startupParam.title = TString("");
   startupParam.lifetimecorrection = 0.0;
-  startupParam.useDKS = false;
+  startupParam.useDKS = MFT_UNDEF;
   startupParam.timeout = 3600;
 }
 
@@ -534,9 +538,9 @@ Int_t musrFT_parse_options(Int_t argc, Char_t *argv[], musrFT_startup_param &sta
         return 2;
       }
       if (!tt.CompareTo("yes", TString::kIgnoreCase))
-        startupParam.useDKS = true;
+        startupParam.useDKS = MFT_DKS;
       else if (!tt.CompareTo("no", TString::kIgnoreCase))
-        startupParam.useDKS = false;
+        startupParam.useDKS = MFT_FFTW;
     } else if (tstr.BeginsWith("--timeout")) {
       if (i+1 >= argc) { // something is wrong since there needs to be an argument here
         cerr << endl << ">> musrFT **ERROR** found option --timeout without argument!" << endl;
@@ -1076,6 +1080,14 @@ Int_t main(Int_t argc, Char_t *argv[])
   }
   startupHandler->SetStartupOptions(startup_options);
 
+  // filter out the use_dks flag
+  if (startupParam.useDKS == MFT_UNDEF) {
+    if (startupHandler->GetStartupOptions()->useDKS)
+      startupParam.useDKS = MFT_DKS;
+    else
+      startupParam.useDKS = MFT_FFTW;
+  }
+
   // defines the raw time-domain data vector
   PPrepFourier data(startupParam.packing, startupParam.bkg_range, startupParam.bkg);
 
@@ -1417,8 +1429,15 @@ Int_t main(Int_t argc, Char_t *argv[])
 
   vector<PFourier*> fourier;
   fourier.resize(histo.size());
+  Bool_t fftw_flag = true;
+  if (startupParam.useDKS == MFT_DKS)
+    fftw_flag = false;
+  if (fftw_flag)
+    cout << "info> will use FFTW" << endl;
+  else
+    cout << "info> will use DKS" << endl;
   for (UInt_t i=0; i<fourier.size(); i++) {
-    fourier[i] = new PFourier(histo[i], unitTag, 0.0, 0.0, true, startupParam.fourierPower, !startupParam.useDKS);
+    fourier[i] = new PFourier(histo[i], unitTag, 0.0, 0.0, true, startupParam.fourierPower, fftw_flag);
   }
 
   // Fourier transform data