first working FFT including DKS.

src/classes/PFourier.cpp

@@ -99,7 +99,6 @@ 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> fStartTime=" << fStartTime << ", fEndTime=" << fEndTime << ", fTimeResolution=" << fTimeResolution << ", fNoOfData=" << fNoOfData << endl;
 
   // check if zero padding is whished
   if (fZeroPaddingPower > 0) {
@@ -247,13 +246,10 @@ void PFourier::Transform(UInt_t apodizationTag)
     int status=0, size=fNoOfBins/2+1;
     // write data to the accelerator
     status=fDks.writeData<double>(fReal_ptr, fInDKS, fNoOfBins);
-    cout << "debug> status=" << status << ": write" << endl;
     // execute the FFT
     status = fDks.callR2CFFT(fReal_ptr, fComp_ptr, 1, dimsize);
-    cout << "debug> status=" << status << ": FFT" << endl;
     // read data from accelerator
     status = fDks.readData< complex<double> >(fComp_ptr, fOutDKS, size);
-    cout << "debug> status=" << status << ": read" << endl;
 #else
     fValid = false;
     return;

src/classes/PMusrCanvas.cpp

@@ -121,6 +121,8 @@ PMusrCanvas::PMusrCanvas()
   fTimeout = 0;
   fTimeoutTimer = 0;
 
+  fStartWithFourier = false;
+  fUseDKS = false;
   fScaleN0AndBkg = true;
   fValid = false;
   fAveragedView     = false;
@@ -185,8 +187,8 @@ PMusrCanvas::PMusrCanvas()
  */
 PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title,
                          Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh,
-                         const Bool_t batch, const Bool_t fourier) :
+                         const Bool_t batch, const Bool_t fourier, const Bool_t useDKS) :
-                         fStartWithFourier(fourier), fBatchMode(batch), fPlotNumber(number)
+                         fStartWithFourier(fourier), fUseDKS(useDKS), fBatchMode(batch), fPlotNumber(number)
 {
   fTimeout = 0;
   fTimeoutTimer = 0;
@@ -240,8 +242,8 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title,
                          Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh,
                          PMsrFourierStructure fourierDefault,
                          const PIntVector markerList, const PIntVector colorList,
-                         const Bool_t batch, const Bool_t fourier) :
+                         const Bool_t batch, const Bool_t fourier, const Bool_t useDKS) :
-                         fStartWithFourier(fourier), fBatchMode(batch),
+                         fStartWithFourier(fourier), fUseDKS(useDKS), fBatchMode(batch),
                          fPlotNumber(number), fFourier(fourierDefault),
                          fMarkerList(markerList), fColorList(colorList)
 {
@@ -3456,7 +3458,12 @@ void PMusrCanvas::HandleFourier()
 
       // calculate fourier transform of the theory
       Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3;
-      PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad);
+      Bool_t useFFTW = true;
+#ifdef HAVE_DKS
+      if ((powerPad >= 20) && fUseDKS)
+        useFFTW = false; // i.e. use DKS
+#endif
+      PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad, useFFTW);
       if (!fourierTheory.IsValid()) {
         cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier theory ..." << endl;
         return;

src/classes/PStartupHandler.cpp

@@ -167,6 +167,7 @@ void PStartupHandler::OnStartDocument()
   fStartupOptions.writeExpectedChisq = false;
   fStartupOptions.estimateN0 = true;
   fStartupOptions.alphaEstimateN0 = 0.0;
+  fStartupOptions.useDKS = false;
 }
 
 //--------------------------------------------------------------------------
@@ -201,6 +202,8 @@ void PStartupHandler::OnStartElement(const Char_t *str, const TList *attributes)
     fKey = eEstimateN0;
   } else if (!strcmp(str, "alpha_estimate_n0")) {
     fKey = eAlphaEstimateN0;
+  } else if (!strcmp(str, "use_dks")) {
+    fKey = eUseDKS;
   } else if (!strcmp(str, "marker")) {
     fKey = eMarker;
   } else if (!strcmp(str, "color")) {
@@ -270,6 +273,11 @@ void PStartupHandler::OnCharacters(const Char_t *str)
       if (tstr.IsFloat())
         fStartupOptions.alphaEstimateN0 = tstr.Atof();
       break;
+    case eUseDKS:
+      tstr = TString(str);
+      if (tstr.BeginsWith("y") || tstr.BeginsWith("Y"))
+        fStartupOptions.useDKS = true;
+      break;
     case eMarker:
       // check that str is a number
       tstr = TString(str);

src/include/PMusr.h

@@ -798,6 +798,7 @@ typedef struct {
   Bool_t writeExpectedChisq; ///< if set to true, expected chisq per block will be written
   Bool_t estimateN0;         ///< if set to true, for single histogram fits N0 will be estimated
   Double_t alphaEstimateN0;  ///< relates the Bkg to N0, i.e. Bkg = alpha*N0
+  Bool_t useDKS;             ///< if set to true, use DKS if present and "sensible"
 } PStartupOptions;
 
 #endif // _PMUSR_H_

src/include/PMusrCanvas.h

@@ -203,12 +203,12 @@ class PMusrCanvas : public TObject, public TQObject
     PMusrCanvas();
     PMusrCanvas(const Int_t number, const Char_t* title,
                 Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, const Bool_t batch,
-                const Bool_t fourier=false);
+                const Bool_t fourier=false, const Bool_t useDKS=false);
     PMusrCanvas(const Int_t number, const Char_t* title,
                 Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh,
                 PMsrFourierStructure fourierDefault,
                 const PIntVector markerList, const PIntVector colorList, const Bool_t batch,
-                const Bool_t fourier=false);
+                const Bool_t fourier=false, const Bool_t useDKS=false);
     virtual ~PMusrCanvas();
 
     virtual Bool_t IsValid() { return fValid; }
@@ -233,6 +233,7 @@ class PMusrCanvas : public TObject, public TQObject
 
   private:
     Bool_t fStartWithFourier; ///< flag if true, the Fourier transform will be presented bypassing the time domain representation
+    Bool_t fUseDKS;           ///< flag if true, use DKS if it is enabled
     Int_t  fTimeout;          ///< timeout after which the Done signal should be emited. If timeout <= 0, no timeout is taking place
     Bool_t fScaleN0AndBkg;    ///< true=N0 and background is scaled to (1/ns), otherwise (1/bin) for the single histogram case
     Bool_t fBatchMode;        ///< musrview in ROOT batch mode

src/include/PStartupHandler.h

@@ -82,7 +82,7 @@ class PStartupHandler : public TObject, public TQObject
     virtual void SetStartupOptions(const PStartupOptions opt) { fStartupOptions = opt; }
 
   private:
-    enum EKeyWords {eEmpty, eComment, eDataPath, eOptions, eWritePerRunBlockChisq, eEstimateN0, eAlphaEstimateN0,
+    enum EKeyWords {eEmpty, eComment, eDataPath, eOptions, eWritePerRunBlockChisq, eEstimateN0, eAlphaEstimateN0, eUseDKS,
                     eFourierSettings, eUnits, eFourierPower, eApodization, ePlot, ePhase, ePhaseIncrement,
                     eRootSettings, eMarkerList, eMarker, 
                     eColorList, eColor};

src/musrFT.cpp

@@ -81,6 +81,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 timeout;                  ///< timeout in (sec) after which musrFT will terminate. if <= 0, no automatic termination will take place.
 } musrFT_startup_param;
 
@@ -148,6 +149,7 @@ void musrFT_syntax()
   cout << endl << "                 provided. This will help on the way to a full fitting model.";
   cout << endl << "    -lc, --lifetimecorrection <fudge>: try to eliminate muon life time decay. Only makes sense for low";
   cout << endl << "                 transverse fields. <fudge> is a tweaking factor and should be kept around 1.0.";
+  cout << endl << "    --useDKS <flag> : if <flag> is true, DKS will be used, otherwise FFTW. Default: false, i.e. FFTW will be used.";
   cout << endl << "    --timeout <timeout> : <timeout> given in seconds after which musrFT terminates.";
   cout << endl << "                 If <timeout> <= 0, no timeout will take place. Default <timeout> is 3600.";
   cout << endl << endl;
@@ -180,6 +182,7 @@ void musrFT_init(musrFT_startup_param &startupParam)
   startupParam.packing = 1;
   startupParam.title = TString("");
   startupParam.lifetimecorrection = 0.0;
+  startupParam.useDKS = false;
   startupParam.timeout = 3600;
 }
 
@@ -519,6 +522,21 @@ Int_t musrFT_parse_options(Int_t argc, Char_t *argv[], musrFT_startup_param &sta
         return 2;
       }
       startupParam.lifetimecorrection = fudge.Atof();
+    } else if (tstr.BeginsWith("--useDKS")) {
+      if (i+1 >= argc) { // something is wrong since there needs to be an argument here
+        cerr << endl << ">> musrFT **ERROR** found option --useDKS without argument!" << endl;
+        return 2;
+      }
+      ++i;
+      TString tt(argv[i]);
+      if (tt.CompareTo("yes", TString::kIgnoreCase) && tt.CompareTo("no", TString::kIgnoreCase)) {
+        cerr << endl << ">> musrFT **ERROR** found option --useDKS with a <flag> which is not yes/no '" << tt << "'." << endl;
+        return 2;
+      }
+      if (!tt.CompareTo("yes", TString::kIgnoreCase))
+        startupParam.useDKS = true;
+      else if (!tt.CompareTo("no", TString::kIgnoreCase))
+        startupParam.useDKS = false;
     } 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;
@@ -712,7 +730,7 @@ Int_t musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, Double_t star
   musrFT_cleanup(hRe);
   for (UInt_t i=1; i<fourierData.size(); i++) {
     hRe = fourierData[i]->GetRealFourier();
-    if (hRe->GetNbinsX()-1 < minSize)
+    if (hRe->GetNbinsX()-1 < (Int_t)minSize)
       minSize = hRe->GetNbinsX()-1;
     musrFT_cleanup(hRe);
   }
@@ -720,7 +738,7 @@ Int_t musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, Double_t star
   for (UInt_t i=0; i<fourierData.size(); i++) {
     hRe = fourierData[i]->GetRealFourier();
     hIm = fourierData[i]->GetImaginaryFourier();
-    for (Int_t j=1; j<minSize; j++) {
+    for (Int_t j=1; j<(Int_t)minSize; j++) {
       dval = hRe->GetBinCenter(j);
       if ((dval >= start) && (dval <= end)) {
         freq.push_back(dval);
@@ -1400,7 +1418,7 @@ Int_t main(Int_t argc, Char_t *argv[])
   vector<PFourier*> fourier;
   fourier.resize(histo.size());
   for (UInt_t i=0; i<fourier.size(); i++) {
-    fourier[i] = new PFourier(histo[i], unitTag, 0.0, 0.0, true, startupParam.fourierPower);
+    fourier[i] = new PFourier(histo[i], unitTag, 0.0, 0.0, true, startupParam.fourierPower, !startupParam.useDKS);
   }
 
   // Fourier transform data

src/musrfit_startup.xml

@@ -18,6 +18,7 @@
       <write_per_run_block_chisq>n</write_per_run_block_chisq>
       <estimate_n0>yes</estimate_n0>
       <alpha_estimate_n0>0.0</alpha_estimate_n0>
+      <use_dks>yes</use_dks>      
     </options>
     <fourier_settings>
         <units>Gauss</units>

src/musrview.cpp

@@ -308,12 +308,12 @@ int main(int argc, char *argv[])
                                      startupHandler->GetMarkerList(),
                                      startupHandler->GetColorList(),
                                      graphicsOutput||asciiOutput,
-                                     fourier);
+                                     fourier, startupHandler->GetStartupOptions()->useDKS);
       else
         musrCanvas = new PMusrCanvas(i, msrHandler->GetMsrTitle()->Data(), 
                                      10+i*100, 10+i*100, 800, 600,
                                      graphicsOutput||asciiOutput,
-                                     fourier);
+                                     fourier, startupHandler->GetStartupOptions()->useDKS);
 
       if (!musrCanvas->IsValid()) {
         cerr << endl << ">> musrview **SEVERE ERROR** Couldn't invoke all necessary objects, will quit.";

src/tests/DKS_FourierTest/dks_fourierTest.cpp

@@ -24,11 +24,13 @@ double millitime()
 //---------------------------------------------------
 void dks_fourierTest_syntax()
 {
-  cout << endl << "usage: dks_fourierTest [useFFTW, N, L | --help] ";
+  cout << endl << "usage: dks_fourierTest [useFFTW N L [dumpAll] | --help] ";
   cout << endl << "       useFFTW : flag, if true -> FFTW, otherwise -> DKS";
   cout << endl << "       N : number of histos";
   cout << endl << "       L : histo length";
-  cout << endl << "       if not given, useFFTW=true, N=8, L=2^20=1048576";
+  cout << endl << "       dumpAll : flag, if true writes the data into a ROOT dump file 'dks_FourierTest.root'";
+  cout << endl << "                 if false, writes only the first and last data set.";
+  cout << endl << "       if not given, useFFTW=true, N=8, L=2^20=1048576, dumpAll=false.";
   cout << endl << "       --help: this help";
   cout << endl << endl;
 }
@@ -42,6 +44,7 @@ int main(int argc, char *argv[])
 
   int N=8, L=1048576;
   Bool_t useFFTW = true;
+  Bool_t dumpAll = false;
 
   switch (argc) {
   case 1:
@@ -58,6 +61,26 @@ int main(int argc, char *argv[])
     N = strtol(argv[2], 0, 10);
     L = strtol(argv[3], 0, 10);
     break;
+  case 5:
+    if (!strcmp(argv[1], "true") || !strcmp(argv[1], "1"))
+      useFFTW = true;
+    else if (!strcmp(argv[1], "false") || !strcmp(argv[1], "0"))
+      useFFTW = false;
+    else {
+      dks_fourierTest_syntax();
+      return 1;
+    }
+    N = strtol(argv[2], 0, 10);
+    L = strtol(argv[3], 0, 10);
+    if (!strcmp(argv[4], "true") || !strcmp(argv[4], "1"))
+      dumpAll = true;
+    else if (!strcmp(argv[4], "false") || !strcmp(argv[4], "0"))
+      dumpAll = false;
+    else {
+      dks_fourierTest_syntax();
+      return 1;
+    }
+    break;
   default:
     dks_fourierTest_syntax();
     return 1;
@@ -124,9 +147,18 @@ int main(int argc, char *argv[])
     cout << endl << "---" << endl;
   }
   TFile fout("dks_fourierTest.root", "recreate");
-  for (unsigned int i=0; i<fourierPowerHistos.size(); i++) {
+  if (dumpAll) {
-    data[i]->Write();
+    for (unsigned int i=0; i<fourierPowerHistos.size(); i++) {
-    fourierPowerHistos[i]->Write();
+      data[i]->Write();
+      fourierPowerHistos[i]->Write();
+    }
+  } else {
+    data[0]->Write();
+    fourierPowerHistos[0]->Write();
+    if (data.size() > 1) {
+      data[data.size()-1]->Write();
+      fourierPowerHistos[fourierPowerHistos.size()-1]->Write();
+    }
   }
   fout.Close();