diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index 39ad4ece..5402c657 100644 --- a/src/classes/PFourier.cpp +++ b/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(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 >(fComp_ptr, fOutDKS, size); - cout << "debug> status=" << status << ": read" << endl; #else fValid = false; return; diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 760d7b2f..df0e8c3e 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/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) : - fStartWithFourier(fourier), fBatchMode(batch), fPlotNumber(number) + const Bool_t batch, const Bool_t fourier, const Bool_t useDKS) : + 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) : - fStartWithFourier(fourier), fBatchMode(batch), + const Bool_t batch, const Bool_t fourier, const Bool_t useDKS) : + 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; diff --git a/src/classes/PStartupHandler.cpp b/src/classes/PStartupHandler.cpp index 3f9b425e..c3d11a8b 100644 --- a/src/classes/PStartupHandler.cpp +++ b/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); diff --git a/src/include/PMusr.h b/src/include/PMusr.h index 8a501d59..f203d717 100644 --- a/src/include/PMusr.h +++ b/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_ diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 19241149..2cddc2b4 100644 --- a/src/include/PMusrCanvas.h +++ b/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 diff --git a/src/include/PStartupHandler.h b/src/include/PStartupHandler.h index f91dc679..71b7c1e4 100644 --- a/src/include/PStartupHandler.h +++ b/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}; diff --git a/src/musrFT.cpp b/src/musrFT.cpp index 37c931c3..e5162e39 100644 --- a/src/musrFT.cpp +++ b/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 : try to eliminate muon life time decay. Only makes sense for low"; cout << endl << " transverse fields. is a tweaking factor and should be kept around 1.0."; + cout << endl << " --useDKS : if is true, DKS will be used, otherwise FFTW. Default: false, i.e. FFTW will be used."; cout << endl << " --timeout : given in seconds after which musrFT terminates."; cout << endl << " If <= 0, no timeout will take place. Default 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 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 &fourierData, Double_t star musrFT_cleanup(hRe); for (UInt_t i=1; iGetRealFourier(); - 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 &fourierData, Double_t star for (UInt_t i=0; iGetRealFourier(); hIm = fourierData[i]->GetImaginaryFourier(); - for (Int_t j=1; jGetBinCenter(j); if ((dval >= start) && (dval <= end)) { freq.push_back(dval); @@ -1400,7 +1418,7 @@ Int_t main(Int_t argc, Char_t *argv[]) vector fourier; fourier.resize(histo.size()); for (UInt_t i=0; i n yes - 0.0 + 0.0 + yes Gauss diff --git a/src/musrview.cpp b/src/musrview.cpp index f91dd6fd..3cdd12c2 100644 --- a/src/musrview.cpp +++ b/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."; diff --git a/src/tests/DKS_FourierTest/dks_fourierTest.cpp b/src/tests/DKS_FourierTest/dks_fourierTest.cpp index 63fe07b5..8e2290b2 100644 --- a/src/tests/DKS_FourierTest/dks_fourierTest.cpp +++ b/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; iWrite(); - fourierPowerHistos[i]->Write(); + if (dumpAll) { + for (unsigned int i=0; iWrite(); + 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();