From 4c24a1e48186c91bf5fdc622536f37cf85694ef6 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 28 Nov 2016 10:13:38 +0100 Subject: [PATCH 01/13] prepare the FOURIER block command 'range_for_phase_correction' for the new FT Re-part automatic phase correction algorithm. --- src/classes/PMsrHandler.cpp | 62 +++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 17 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 9bd06d26..f5b75804 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -3844,7 +3844,7 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) TObjArray *tokens = 0; TObjString *ostr = 0; - TString str; + TString str, pcStr=TString(""); Int_t ival; @@ -3993,21 +3993,10 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) } } } else if (iter->fLine.BeginsWith("range_for_phase_correction", TString::kIgnoreCase)) { - if (tokens->GetEntries() < 3) { // range values are missing - error = true; - continue; - } else { - for (UInt_t i=0; i<2; i++) { - ostr = dynamic_cast(tokens->At(i+1)); - str = ostr->GetString(); - if (str.IsFloat()) { - fourier.fRangeForPhaseCorrection[i] = str.Atof(); - } else { - error = true; - continue; - } - } - } + // keep the string. It can only be handled at the very end of the FOURIER block evaluation + // since it needs potentially the range input and there is no guaranty this is already + // available at this point. + pcStr = iter->fLine; } else if (iter->fLine.BeginsWith("range", TString::kIgnoreCase)) { // fourier plot range if (tokens->GetEntries() < 3) { // plot range values are missing error = true; @@ -4045,6 +4034,45 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) tokens = 0; } + // handle range_for_phase_correction if present + if ((pcStr.Length() != 0) && !error) { + // tokenize line + tokens = pcStr.Tokenize(" \t"); + + switch (tokens->GetEntries()) { + case 2: + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (!str.CompareTo("all", TString::kIgnoreCase)) { + fourier.fRangeForPhaseCorrection[0] = fourier.fPlotRange[0]; + fourier.fRangeForPhaseCorrection[1] = fourier.fPlotRange[1]; + } else { + error = true; + } + break; + case 3: + for (UInt_t i=0; i<2; i++) { + ostr = dynamic_cast(tokens->At(i+1)); + str = ostr->GetString(); + if (str.IsFloat()) { + fourier.fRangeForPhaseCorrection[i] = str.Atof(); + } else { + error = true; + } + } + break; + default: + error = true; + break; + } + + // clean up + if (tokens) { + delete tokens; + tokens = 0; + } + } + if (error) { cerr << endl << ">> PMsrHandler::HandleFourierEntry: **ERROR** in line " << iter->fLineNo << ":"; cerr << endl; @@ -4061,7 +4089,7 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) cerr << endl << ">> [apodization none | weak | medium | strong]"; cerr << endl << ">> [plot real | imag | real_and_imag | power | phase]"; cerr << endl << ">> [phase value]"; - cerr << endl << ">> [range_for_phase_correction min max]"; + cerr << endl << ">> [range_for_phase_correction min max | all]"; cerr << endl << ">> [range min max]"; cerr << endl; } else { // save last run found From 1e5bfb821682c617706f1fe848cf16170aac3c2e Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 28 Nov 2016 11:02:03 +0100 Subject: [PATCH 02/13] when activating the cross hair cursor in the musrview canvas, enable the event status bar as well. This allows the user to see exactly where the cursor is. This is for instance helpful for Fourier spectra. --- src/classes/PMusrCanvas.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 9a0a3b49..44d31939 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1142,10 +1142,13 @@ void PMusrCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected) } } else if (x == 'c') { Int_t state = fDataTheoryPad->GetCrosshair(); - if (state == 0) + if (state == 0) { + fMainCanvas->ToggleEventStatus(); fDataTheoryPad->SetCrosshair(2); - else + } else { + fMainCanvas->ToggleEventStatus(); fDataTheoryPad->SetCrosshair(0); + } fMainCanvas->Update(); } else { fMainCanvas->Update(); From f7834aaead4696e2e295c7fb821bff66e8b75e14 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Tue, 29 Nov 2016 11:16:45 +0100 Subject: [PATCH 03/13] added class PFTPhaseCorrection which allows to calculate the phase optimzied real part of the Fourier transform. musrFT is already using it, for musrview it is still missing. --- src/classes/PFourier.cpp | 287 +++++++++++++++++++++++++++++---- src/classes/PFourierCanvas.cpp | 55 ++++--- src/include/PFourier.h | 55 ++++++- src/include/PFourierCanvas.h | 2 +- 4 files changed, 339 insertions(+), 60 deletions(-) diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index ec516096..b6034b26 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -37,7 +37,9 @@ using namespace std; #include "TF1.h" #include "TAxis.h" -//#include "TFile.h" +#include "Minuit2/FunctionMinimum.h" +#include "Minuit2/MnUserParameters.h" +#include "Minuit2/MnMinimize.h" #include "PMusr.h" #include "PFourier.h" @@ -45,6 +47,227 @@ using namespace std; #define PI 3.14159265358979312 #define PI_HALF 1.57079632679489656 +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// PFTPhaseCorrection +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + * + */ +PFTPhaseCorrection::PFTPhaseCorrection(const Int_t minBin, const Int_t maxBin) : + fMinBin(minBin), fMaxBin(maxBin) +{ + Init(); +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + * + */ +PFTPhaseCorrection::PFTPhaseCorrection(vector &reFT, vector &imFT, const Int_t minBin, const Int_t maxBin) : + fReal(reFT), fImag(imFT), fMinBin(minBin), fMaxBin(maxBin) +{ + Init(); + + if (fMinBin == -1) + fMinBin = 0; + if (fMaxBin == -1) + fMaxBin = fReal.size(); + if (fMaxBin > fReal.size()) + fMaxBin = fReal.size(); + + fRealPh.resize(fReal.size()); + fImagPh.resize(fReal.size()); +} + +//-------------------------------------------------------------------------- +// Minimize (public) +//-------------------------------------------------------------------------- +/** + * + */ +void PFTPhaseCorrection::Minimize() +{ + // create Minuit2 parameters + ROOT::Minuit2::MnUserParameters upar; + + upar.Add("c0", fPh_c0, 0.5); + upar.Add("c1", fPh_c1, 0.5); + + // create minimizer + ROOT::Minuit2::MnMinimize mn_min(*this, upar); + + // minimize + ROOT::Minuit2::FunctionMinimum min = mn_min(); + + if (min.IsValid()) { + fPh_c0 = min.UserState().Value("c0"); + fPh_c1 = min.UserState().Value("c1"); + fMin = min.Fval(); + } else { + fMin = -1.0; + fValid = false; + cout << endl << ">> **WARNING** minimize failed to find a minimum for the real FT phase correction ..." << endl; + return; + } +} + +//-------------------------------------------------------------------------- +// GetPhaseCorrectionParam (public) +//-------------------------------------------------------------------------- +/** + * + */ +Double_t PFTPhaseCorrection::GetPhaseCorrectionParam(UInt_t idx) +{ + Double_t result=0.0; + + if (idx == 0) + result = fPh_c0; + else if (idx == 1) + result = fPh_c1; + else + cerr << ">> **ERROR** requested phase correction parameter with index=" << idx << " does not exist!" << endl; + + return result; +} + +//-------------------------------------------------------------------------- +// GetMinimum (public) +//-------------------------------------------------------------------------- +/** + * + */ +Double_t PFTPhaseCorrection::GetMinimum() +{ + if (!fValid) { + cerr << ">> **ERROR** requested minimum is invalid!" << endl; + return -1.0; + } + + return fMin; +} + +//-------------------------------------------------------------------------- +// Init (private) +//-------------------------------------------------------------------------- +/** + * + */ +void PFTPhaseCorrection::Init() +{ + fValid = true; + fPh_c0 = 0.0; + fPh_c1 = 0.0; + fGamma = 1.0; + fMin = -1.0; +} + +//-------------------------------------------------------------------------- +// CalcPhasedFT (private) +//-------------------------------------------------------------------------- +/** + * + */ +void PFTPhaseCorrection::CalcPhasedFT() const +{ + Double_t phi=0.0; + Double_t w=0.0; + + for (UInt_t i=0; i 1.0e-15) { + hh = dval / norm; + entropy -= hh * log(hh); + } + } + + return entropy; +} + +//-------------------------------------------------------------------------- +// operator() (private) +//-------------------------------------------------------------------------- +/** + * + */ +double PFTPhaseCorrection::operator()(const vector &par) const +{ + // par : [0]: c0, [1]: c1 + + fPh_c0 = par[0]; + fPh_c1 = par[1]; + + CalcPhasedFT(); + CalcRealPhFTDerivative(); + + return Entropy()+Penalty(); +} + +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// PFourier //-------------------------------------------------------------------------- // Constructor //-------------------------------------------------------------------------- @@ -73,6 +296,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi fValid = true; fIn = 0; fOut = 0; + fPhCorrectedReFT = 0; fApodization = F_APODIZATION_NONE; @@ -160,6 +384,8 @@ PFourier::~PFourier() fftw_free(fIn); if (fOut) fftw_free(fOut); + if (fPhCorrectedReFT) + delete fPhCorrectedReFT; } //-------------------------------------------------------------------------- @@ -262,26 +488,19 @@ TH1F* PFourier::GetRealFourier(const Double_t scale) // GetPhaseOptRealFourier (public) //-------------------------------------------------------------------------- /** - *

returns the phase corrected real Fourier transform. The correction angle is - * past back as well. - *

Currently it simply does the following thing: (i) rotate the complex Fourier - * transform through all angles in 1/2° steps, i.e. - * \f$ f_{\rm rot} = (f_{\rm real} + i f_{\rm imag}) \exp(- \alpha)\f$, - * hence \f$ f_{\rm rot} = f_{\rm real} \cos(\alpha) + f_{\rm imag} \sin(\alpha)\f$. - * (ii) search the maximum of \f$ sum_{\alpha} {\cal R}\{f_{\rm rot}\}\f$ for all - * \f$\alpha\f$. From this one gets \f$\alpha_{\rm opt}\f$. (iii) The 'optimal' - * real Fourier transform is \f$f_{\rm opt} = (f_{\rm real} + i f_{\rm imag}) - * \exp(- \alpha_{\rm opt})\f$. + *

returns the phase corrected real Fourier transform. * * \return the TH1F histo of the phase 'optimzed' real Fourier transform. * - * \param phase return value of the optimal phase + * \param phase return value of the optimal phase dispersion phase[0]+phase[1]*i/N * \param scale normalisation factor * \param min minimal freq / field from which to optimise. Given in the choosen unit. * \param max maximal freq / field up to which to optimise. Given in the choosen unit. */ -TH1F* PFourier::GetPhaseOptRealFourier(Double_t &phase, const Double_t scale, const Double_t min, const Double_t max) +TH1F* PFourier::GetPhaseOptRealFourier(vector &phase, const Double_t scale, const Double_t min, const Double_t max) { + phase.resize(2); // c0, c1 + UInt_t noOfFourierBins = 0; if (fNoOfBins % 2 == 0) noOfFourierBins = fNoOfBins/2; @@ -312,32 +531,28 @@ TH1F* PFourier::GetPhaseOptRealFourier(Double_t &phase, const Double_t scale, co imagF.push_back(fOut[i][1]); } - // rotate trough real/imag Fourier through 360° with a 1/2° resolution and keep the integral - Double_t rotRealIntegral[720]; - Double_t sum = 0.0; - Double_t da = 8.72664625997164774e-03; // pi / 360 - for (UInt_t i=0; i<720; i++) { - sum = 0.0; - for (UInt_t j=0; jMinimize(); + if (!fPhCorrectedReFT->IsValid()) { + fValid = false; + cerr << endl << "**ERROR** could not fina a valid phase correction minimum." << endl; + return 0; } - - // keep the optimal phase - phase = maxRotBin*da; + phase[0] = fPhCorrectedReFT->GetPhaseCorrectionParam(0); + phase[1] = fPhCorrectedReFT->GetPhaseCorrectionParam(1); // clean up + if (fPhCorrectedReFT) { + delete fPhCorrectedReFT; + fPhCorrectedReFT = 0; + } realF.clear(); imagF.clear(); @@ -355,8 +570,10 @@ TH1F* PFourier::GetPhaseOptRealFourier(Double_t &phase, const Double_t scale, co } // fill realFourier vector + Double_t ph; for (UInt_t i=0; iSetBinContent(i+1, scale*(fOut[i][0]*cos(phase) + fOut[i][1]*sin(phase))); + ph = phase[0] + phase[1] * (Double_t)((Int_t)i-(Int_t)minBin) / (Double_t)(maxBin-minBin); + realPhaseOptFourier->SetBinContent(i+1, scale*(fOut[i][0]*cos(ph) - fOut[i][1]*sin(ph))); realPhaseOptFourier->SetBinError(i+1, 0.0); } diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp index 07813ff4..0b9cbc1f 100644 --- a/src/classes/PFourierCanvas.cpp +++ b/src/classes/PFourierCanvas.cpp @@ -624,6 +624,12 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName) xMinBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetFirst(); xMaxBin = fFourierHistos[0].dataFourierIm->GetXaxis()->GetLast(); break; + case FOURIER_PLOT_REAL_AND_IMAG: + xAxis = fFourierHistos[0].dataFourierRe->GetXaxis()->GetTitle(); + yAxis = TString("Real+Imag"); + xMinBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetFirst(); + xMaxBin = fFourierHistos[0].dataFourierRe->GetXaxis()->GetLast(); + break; case FOURIER_PLOT_POWER: xAxis = fFourierHistos[0].dataFourierPwr->GetXaxis()->GetTitle(); yAxis = TString("Power"); @@ -703,42 +709,46 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName) for (UInt_t j=0; jGetBinCenter(i) << ", " << fFourierHistos[j].dataFourierRe->GetBinContent(i) << ", "; - break; + fout << fFourierHistos[j].dataFourierRe->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierRe->GetBinContent(i) << ", "; + break; case FOURIER_PLOT_IMAG: - fout << fFourierHistos[j].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierIm->GetBinContent(i) << ", "; - break; + fout << fFourierHistos[j].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierIm->GetBinContent(i) << ", "; + break; + case FOURIER_PLOT_REAL_AND_IMAG: + break; case FOURIER_PLOT_POWER: - fout << fFourierHistos[j].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPwr->GetBinContent(i) << ", "; - break; + fout << fFourierHistos[j].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPwr->GetBinContent(i) << ", "; + break; case FOURIER_PLOT_PHASE: - fout << fFourierHistos[j].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhase->GetBinContent(i) << ", "; - break; + fout << fFourierHistos[j].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhase->GetBinContent(i) << ", "; + break; case FOURIER_PLOT_PHASE_OPT_REAL: - fout << fFourierHistos[j].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhaseOptReal->GetBinContent(i) << ", "; - break; + fout << fFourierHistos[j].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[j].dataFourierPhaseOptReal->GetBinContent(i) << ", "; + break; default: - break; + break; } } switch (fCurrentPlotView) { case FOURIER_PLOT_REAL: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinContent(i) << endl; - break; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierRe->GetBinContent(i) << endl; + break; case FOURIER_PLOT_IMAG: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinContent(i) << endl; - break; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierIm->GetBinContent(i) << endl; + break; + case FOURIER_PLOT_REAL_AND_IMAG: + break; case FOURIER_PLOT_POWER: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinContent(i) << endl; - break; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPwr->GetBinContent(i) << endl; + break; case FOURIER_PLOT_PHASE: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinContent(i) << endl; - break; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhase->GetBinContent(i) << endl; + break; case FOURIER_PLOT_PHASE_OPT_REAL: - fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinContent(i) << endl; - break; + fout << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierHistos[fFourierHistos.size()-1].dataFourierPhaseOptReal->GetBinContent(i) << endl; + break; default: - break; + break; } } } @@ -805,7 +815,6 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierPwr = fFourier[i]->GetPowerFourier(); fFourierHistos[i].dataFourierPhase = fFourier[i]->GetPhaseFourier(); fFourierHistos[i].dataFourierPhaseOptReal = fFourier[i]->GetPhaseOptRealFourier(fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); - cout << "debug> histo[" << i << "]: opt. phase = " << fFourierHistos[i].optPhase * 180.0 / TMath::Pi() << "°" << endl; } // rescale histo to abs(maximum) == 1 diff --git a/src/include/PFourier.h b/src/include/PFourier.h index bff94450..38f7e7ee 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -30,8 +30,13 @@ #ifndef _PFOURIER_H_ #define _PFOURIER_H_ +#include +using namespace std; + #include "fftw3.h" +#include "Minuit2/FCNBase.h" + #include "PMusr.h" #define F_APODIZATION_NONE 1 @@ -39,6 +44,52 @@ #define F_APODIZATION_MEDIUM 3 #define F_APODIZATION_STRONG 4 +/** + * Re Fourier phase correction + */ +class PFTPhaseCorrection : public ROOT::Minuit2::FCNBase +{ + public: + PFTPhaseCorrection(const Int_t minBin=-1, const Int_t maxBin=-1); + PFTPhaseCorrection(vector &reFT, vector &imFT, const Int_t minBin=-1, const Int_t maxBin=-1); + virtual ~PFTPhaseCorrection() {} + + virtual Bool_t IsValid() { return fValid; } + virtual void Minimize(); + + virtual void SetGamma(const Double_t gamma) { fGamma = gamma; } + virtual void SetPh(const Double_t c0, const Double_t c1) { fPh_c0 = c0; fPh_c1 = c1; CalcPhasedFT(); CalcRealPhFTDerivative(); } + + virtual Double_t GetGamma() { return fGamma; } + virtual Double_t GetPhaseCorrectionParam(UInt_t idx); + virtual Double_t GetMinimum(); + + private: + Bool_t fValid; + + vector fReal; /// original real Fourier data set + vector fImag; /// original imag Fourier data set + mutable vector fRealPh; /// phased real Fourier data set + mutable vector fImagPh; /// phased imag Fourier data set + mutable vector fRealPhD; /// 1st derivative of fRealPh + + Double_t fMinBin; /// minimum bin from Fourier range to be used for the phase correction estimate + Double_t fMaxBin; /// maximum bin from Fourier range to be used for the phase correction estimate + mutable Double_t fPh_c0; /// constant part of the phase dispersion used for the phase correction + mutable Double_t fPh_c1; /// linear part of the phase dispersion used for the phase correction + Double_t fGamma; /// gamma parameter to balance between entropy and penalty + Double_t fMin; /// keeps the minimum of the entropy/penalty minimization + + virtual void Init(); + virtual void CalcPhasedFT() const; + virtual void CalcRealPhFTDerivative() const; + virtual Double_t Penalty() const; + virtual Double_t Entropy() const; + + virtual Double_t Up() const { return 1.0; } + virtual Double_t operator()(const vector&) const; +}; + /** * muSR Fourier class. */ @@ -57,7 +108,7 @@ class PFourier virtual Double_t GetResolution() { return fResolution; } virtual Double_t GetMaxFreq(); virtual TH1F* GetRealFourier(const Double_t scale = 1.0); - virtual TH1F* GetPhaseOptRealFourier(Double_t &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); + virtual TH1F* GetPhaseOptRealFourier(vector &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0); virtual TH1F* GetPowerFourier(const Double_t scale = 1.0); virtual TH1F* GetPhaseFourier(const Double_t scale = 1.0); @@ -85,6 +136,8 @@ class PFourier fftw_complex *fIn; ///< real part of the Fourier transform fftw_complex *fOut; ///< imaginary part of the Fourier transform + PFTPhaseCorrection *fPhCorrectedReFT; + virtual void PrepareFFTwInputData(UInt_t apodizationTag); virtual void ApodizeData(Int_t apodizationTag); }; diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h index 1229b34c..3aff49e9 100644 --- a/src/include/PFourierCanvas.h +++ b/src/include/PFourierCanvas.h @@ -71,7 +71,7 @@ typedef struct { TH1F *dataFourierPwr; ///< power spectrum of the Fourier transform of the data histogram TH1F *dataFourierPhase; ///< phase spectrum of the Fourier transform of the data histogram TH1F *dataFourierPhaseOptReal; ///< phase otpimized real Fourier transform of the data histogram - Double_t optPhase; ///< optimal phase which maximizes the real Fourier + vector optPhase; ///< optimal phase which maximizes the real Fourier } PFourierCanvasDataSet; //------------------------------------------------------------------------ From a361c7f7bb44e4deb3ec409d7755eade5702f92b Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Tue, 29 Nov 2016 21:06:10 +0100 Subject: [PATCH 04/13] added locally the LD_LIBRARY_PATH. This is needed for macOS 10.12 if ROOT is installed in a local user directory. --- src/musredit_qt5/PFitOutputHandler.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/musredit_qt5/PFitOutputHandler.cpp b/src/musredit_qt5/PFitOutputHandler.cpp index 28ca0298..277618d6 100644 --- a/src/musredit_qt5/PFitOutputHandler.cpp +++ b/src/musredit_qt5/PFitOutputHandler.cpp @@ -62,6 +62,10 @@ PFitOutputHandler::PFitOutputHandler(QString workingDirectory, QVector // QProcess related code fProc = new QProcess( this ); + // make sure that the system environment variables are properly set + QProcessEnvironment env = QProcessEnvironment::systemEnvironment(); + env.insert("LD_LIBRARY_PATH", env.value("ROOTSYS") + "/lib:" + env.value("LD_LIBRARY_PATH")); + fProc->setProcessEnvironment(env); fProc->setWorkingDirectory(workingDirectory); // Set up the command and arguments. From b106ac61a110d7fb0ab2f498a3e144615c998ec5 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Wed, 7 Dec 2016 17:02:18 +0100 Subject: [PATCH 05/13] added phase optimized Re FT for musrview. Still not perfect. --- src/classes/PMsrHandler.cpp | 12 +- src/classes/PMusrCanvas.cpp | 366 ++++++++++++++++++++++++++++-------- src/include/PMusrCanvas.h | 32 ++-- 3 files changed, 318 insertions(+), 92 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index f5b75804..d126ccd9 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1113,8 +1113,10 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout << "POWER"; } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE) { fout << "PHASE"; + } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE_OPT_REAL) { + fout << "PHASE_OPT_REAL"; } - fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE"; + fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL"; fout << endl; } else if (sstr.BeginsWith("phase")) { if (fFourier.fPhaseParamNo > 0) { @@ -2146,8 +2148,10 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co fout << "POWER"; } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE) { fout << "PHASE"; + } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE_OPT_REAL) { + fout << "PHASE_OPT_REAL"; } - fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE"; + fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL"; fout << endl; } @@ -3951,6 +3955,8 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) fourier.fPlotTag = FOURIER_PLOT_POWER; } else if (!str.CompareTo("phase", TString::kIgnoreCase)) { fourier.fPlotTag = FOURIER_PLOT_PHASE; + } else if (!str.CompareTo("phase_opt_real", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_PHASE_OPT_REAL; } else { // unrecognized plot tag error = true; continue; @@ -4087,7 +4093,7 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) cerr << endl << ">> 0 <= n <= 20 are allowed values"; cerr << endl << ">> [dc-corrected true | false]"; cerr << endl << ">> [apodization none | weak | medium | strong]"; - cerr << endl << ">> [plot real | imag | real_and_imag | power | phase]"; + cerr << endl << ">> [plot real | imag | real_and_imag | power | phase | phase_opt_real]"; cerr << endl << ">> [phase value]"; cerr << endl << ">> [range_for_phase_correction min max | all]"; cerr << endl << ">> [range min max]"; diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 44d31939..2cae4d9e 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -815,6 +815,12 @@ void PMusrCanvas::UpdateDataTheoryPad() fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); } break; + case FOURIER_PLOT_PHASE_OPT_REAL: + fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + if (!fBatchMode) { + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL); + } + break; default: fCurrentPlotView = PV_FOURIER_PWR; if (!fBatchMode) { @@ -1096,6 +1102,11 @@ void PMusrCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected) fCurrentPlotView = PV_FOURIER_PHASE; fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); break; + case FOURIER_PLOT_PHASE_OPT_REAL: + fPreviousPlotView = fCurrentPlotView; + fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL); + break; default: break; } @@ -1342,6 +1353,29 @@ void PMusrCanvas::HandleMenuPopup(Int_t id) HandleFourierDifference(); PlotFourierDifference(); } + } else if (id == P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL) { + // set appropriate plot view + fPreviousPlotView = fCurrentPlotView; + fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + // uncheck data + fPopupMain->UnCheckEntry(P_MENU_ID_DATA+P_MENU_PLOT_OFFSET*fPlotNumber); + // check appropriate fourier popup item + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(id); + // enable phase increment/decrement + fPopupFourier->EnableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS); + fPopupFourier->EnableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS); + // handle fourier phase + if (!fDifferenceView) { + HandleFourier(); + PlotFourier(); + } else { + if (previousPlotView == PV_DATA) + HandleDifferenceFourier(); + else + HandleFourierDifference(); + PlotFourierDifference(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS) { IncrementFourierPhase(); } else if (id == P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS) { @@ -1749,6 +1783,24 @@ void PMusrCanvas::ExportData(const Char_t *fileName) } } break; + case PV_FOURIER_PHASE_OPT_REAL: + // get current x-range + xminBin = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierPhaseOptReal, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierPhaseOptReal, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iAddEntry("Show Real+Imag", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_REAL_AND_IMAG); fPopupFourier->AddEntry("Show Power", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PWR); fPopupFourier->AddEntry("Show Phase", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); + fPopupFourier->AddEntry("Show PhaseOptReal", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL); fPopupFourier->AddSeparator(); fPopupFourier->AddEntry("Phase +", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS); fPopupFourier->AddEntry("Phase -", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS); @@ -2417,16 +2473,19 @@ void PMusrCanvas::InitDataSet(PMusrCanvasDataSet &dataSet) dataSet.dataFourierIm = 0; dataSet.dataFourierPwr = 0; dataSet.dataFourierPhase = 0; + dataSet.dataFourierPhaseOptReal = 0; dataSet.theory = 0; dataSet.theoryFourierRe = 0; dataSet.theoryFourierIm = 0; dataSet.theoryFourierPwr = 0; dataSet.theoryFourierPhase = 0; + dataSet.theoryFourierPhaseOptReal = 0; dataSet.diff = 0; dataSet.diffFourierRe = 0; dataSet.diffFourierIm = 0; dataSet.diffFourierPwr = 0; dataSet.diffFourierPhase = 0; + dataSet.diffFourierPhaseOptReal = 0; dataSet.dataRange = 0; } @@ -2488,6 +2547,10 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet) delete dataSet.dataFourierPhase; dataSet.dataFourierPhase = 0; } + if (dataSet.dataFourierPhaseOptReal) { + delete dataSet.dataFourierPhaseOptReal; + dataSet.dataFourierPhaseOptReal = 0; + } if (dataSet.theory) { delete dataSet.theory; dataSet.theory = 0; @@ -2508,6 +2571,10 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet) delete dataSet.theoryFourierPhase; dataSet.theoryFourierPhase = 0; } + if (dataSet.theoryFourierPhaseOptReal) { + delete dataSet.theoryFourierPhaseOptReal; + dataSet.theoryFourierPhaseOptReal = 0; + } if (dataSet.diff) { delete dataSet.diff; dataSet.diff = 0; @@ -2528,6 +2595,10 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet) delete dataSet.diffFourierPhase; dataSet.diffFourierPhase = 0; } + if (dataSet.diffFourierPhaseOptReal) { + delete dataSet.diffFourierPhaseOptReal; + dataSet.diffFourierPhaseOptReal = 0; + } if (dataSet.dataRange) { delete dataSet.dataRange; dataSet.dataRange = 0; @@ -3260,6 +3331,11 @@ void PMusrCanvas::HandleDifference() */ void PMusrCanvas::HandleFourier() { + PDoubleVector phaseParam; + Double_t min=-1.0, max=-1.0; + Double_t re, im, ph; + char hName[1024]; + // check if plot type is appropriate for fourier if (fPlotType == MSR_PLOT_NON_MUSR) return; @@ -3269,7 +3345,7 @@ void PMusrCanvas::HandleFourier() Int_t bin; double startTime = fXmin; double endTime = fXmax; - if (!fStartWithFourier) { // fHistoFrame presen, hence get start/end from it + if (!fStartWithFourier) { // fHistoFrame present, hence get start/end from it bin = fHistoFrame->GetXaxis()->GetFirst(); startTime = fHistoFrame->GetBinCenter(bin); bin = fHistoFrame->GetXaxis()->GetLast(); @@ -3294,6 +3370,12 @@ void PMusrCanvas::HandleFourier() // get phase part of the data fData[i].dataFourierPhase = fourierData.GetPhaseFourier(); + // get phase optimized real fourier from data + min = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; + max = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; + // eventually min/max needs to be extracted from the 'range_for_phase_correction' parameter of the msr-file + fData[i].dataFourierPhaseOptReal = fourierData.GetPhaseOptRealFourier(phaseParam, scale, min, max); + // set marker and line color fData[i].dataFourierRe->SetMarkerColor(fData[i].data->GetMarkerColor()); fData[i].dataFourierRe->SetLineColor(fData[i].data->GetLineColor()); @@ -3303,17 +3385,22 @@ void PMusrCanvas::HandleFourier() fData[i].dataFourierPwr->SetLineColor(fData[i].data->GetLineColor()); fData[i].dataFourierPhase->SetMarkerColor(fData[i].data->GetMarkerColor()); fData[i].dataFourierPhase->SetLineColor(fData[i].data->GetLineColor()); + fData[i].dataFourierPhaseOptReal->SetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPhaseOptReal->SetLineColor(fData[i].data->GetLineColor()); // set marker size fData[i].dataFourierRe->SetMarkerSize(1); fData[i].dataFourierIm->SetMarkerSize(1); fData[i].dataFourierPwr->SetMarkerSize(1); fData[i].dataFourierPhase->SetMarkerSize(1); + fData[i].dataFourierPhaseOptReal->SetMarkerSize(1); + // set marker type fData[i].dataFourierRe->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierIm->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPwr->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPhase->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + fData[i].dataFourierPhaseOptReal->SetMarkerStyle(fData[i].data->GetMarkerStyle()); // calculate fourier transform of the theory Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3; @@ -3333,16 +3420,33 @@ void PMusrCanvas::HandleFourier() // get phase part of the data fData[i].theoryFourierPhase = fourierTheory.GetPhaseFourier(); + // get phase optimized real fourier from data + + // clone theory Re FT + strcpy(hName, fData[i].theoryFourierPhase->GetName()); + strcat(hName, "_Opt_Real"); + fData[i].theoryFourierPhaseOptReal = (TH1F*) fData[i].theoryFourierRe->Clone(hName); + + // rotate the theory according to the optimized phase parameters + // first find minBin for min of the phase correction + Int_t minBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindBin(min); + Int_t maxBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindBin(max); + for (Int_t j=1; jGetNbinsX(); j++) { + ph = phaseParam[0] + phaseParam[1] * (Double_t)(j-minBin+1) / (Double_t)(maxBin-minBin); + re = fData[i].theoryFourierRe->GetBinContent(j) * cos(ph) - fData[i].theoryFourierIm->GetBinContent(j) * sin(ph); + fData[i].theoryFourierPhaseOptReal->SetBinContent(j, re); + } + // set line colors for the theory fData[i].theoryFourierRe->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierIm->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPwr->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPhase->SetLineColor(fData[i].theory->GetLineColor()); + fData[i].theoryFourierPhaseOptReal->SetLineColor(fData[i].theory->GetLineColor()); } // apply global phase if present if (fFourier.fPhase != 0.0) { - double re, im; const double cp = TMath::Cos(fFourier.fPhase/180.0*TMath::Pi()); const double sp = TMath::Sin(fFourier.fPhase/180.0*TMath::Pi()); @@ -3372,6 +3476,7 @@ void PMusrCanvas::HandleFourier() } } +/* as: will be obsolate ... // find optimal Fourier phase if range is given if ((fFourier.fRangeForPhaseCorrection[0] != -1.0) && (fFourier.fRangeForPhaseCorrection[1] != -1.0)) { @@ -3405,6 +3510,7 @@ void PMusrCanvas::HandleFourier() } } } +*/ } } @@ -3546,6 +3652,12 @@ void PMusrCanvas::HandleFourierDifference() fData[i].dataFourierPhase->GetXaxis()->GetXmin(), fData[i].dataFourierPhase->GetXaxis()->GetXmax()); + // phase optimized real part + name = TString(fData[i].dataFourierPhaseOptReal->GetTitle()) + "_diff"; + fData[i].diffFourierPhaseOptReal = new TH1F(name, name, fData[i].dataFourierPhaseOptReal->GetNbinsX(), + fData[i].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[i].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + // calculate difference for (UInt_t j=1; jGetEntries(); j++) { dvalx = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j @@ -3564,6 +3676,10 @@ void PMusrCanvas::HandleFourierDifference() theoBin = fData[i].theoryFourierPhase->FindBin(dvalx); // get the theory x-axis bin dval = fData[i].dataFourierPhase->GetBinContent(j) - fData[i].theoryFourierPhase->GetBinContent(theoBin); fData[i].diffFourierPhase->SetBinContent(j, dval); + dvalx = fData[i].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierPhaseOptReal->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierPhaseOptReal->GetBinContent(j) - fData[i].theoryFourierPhaseOptReal->GetBinContent(theoBin); + fData[i].diffFourierPhaseOptReal->SetBinContent(j, dval); } } @@ -3577,17 +3693,21 @@ void PMusrCanvas::HandleFourierDifference() fData[i].diffFourierPwr->SetLineColor(fData[i].dataFourierPwr->GetLineColor()); fData[i].diffFourierPhase->SetMarkerColor(fData[i].dataFourierPhase->GetMarkerColor()); fData[i].diffFourierPhase->SetLineColor(fData[i].dataFourierPhase->GetLineColor()); + fData[i].diffFourierPhaseOptReal->SetMarkerColor(fData[i].dataFourierPhaseOptReal->GetMarkerColor()); + fData[i].diffFourierPhaseOptReal->SetLineColor(fData[i].dataFourierPhaseOptReal->GetLineColor()); // set marker size fData[i].diffFourierRe->SetMarkerSize(1); fData[i].diffFourierIm->SetMarkerSize(1); fData[i].diffFourierPwr->SetMarkerSize(1); fData[i].diffFourierPhase->SetMarkerSize(1); + fData[i].diffFourierPhaseOptReal->SetMarkerSize(1); // set marker type fData[i].diffFourierRe->SetMarkerStyle(fData[i].dataFourierRe->GetMarkerStyle()); fData[i].diffFourierIm->SetMarkerStyle(fData[i].dataFourierIm->GetMarkerStyle()); fData[i].diffFourierPwr->SetMarkerStyle(fData[i].dataFourierPwr->GetMarkerStyle()); fData[i].diffFourierPhase->SetMarkerStyle(fData[i].dataFourierPhase->GetMarkerStyle()); + fData[i].diffFourierPhaseOptReal->SetMarkerStyle(fData[i].dataFourierPhaseOptReal->GetMarkerStyle()); // set diffFourierTag fData[i].diffFourierTag = 2; // f-d @@ -3646,6 +3766,12 @@ void PMusrCanvas::HandleAverage() fData[0].dataFourierPhase->GetXaxis()->GetXmin(), fData[0].dataFourierPhase->GetXaxis()->GetXmax()); } + if (fData[0].dataFourierPhaseOptReal != 0) { + name = TString(fData[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.dataFourierPhaseOptReal = new TH1F(name, name, fData[0].dataFourierPhaseOptReal->GetNbinsX(), + fData[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } if (fData[0].theory != 0) { name = TString(fData[0].theory->GetTitle()) + "_avg"; fDataAvg.theory = new TH1F(name, name, fData[0].theory->GetNbinsX(), @@ -3676,6 +3802,12 @@ void PMusrCanvas::HandleAverage() fData[0].theoryFourierPhase->GetXaxis()->GetXmin(), fData[0].theoryFourierPhase->GetXaxis()->GetXmax()); } + if (fData[0].theoryFourierPhaseOptReal != 0) { + name = TString(fData[0].theoryFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.theoryFourierPhaseOptReal = new TH1F(name, name, fData[0].theoryFourierPhaseOptReal->GetNbinsX(), + fData[0].theoryFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].theoryFourierPhaseOptReal->GetXaxis()->GetXmax()); + } if (fData[0].diff != 0) { name = TString(fData[0].diff->GetTitle()) + "_avg"; fDataAvg.diff = new TH1F(name, name, fData[0].diff->GetNbinsX(), @@ -3706,6 +3838,12 @@ void PMusrCanvas::HandleAverage() fData[0].diffFourierPhase->GetXaxis()->GetXmin(), fData[0].diffFourierPhase->GetXaxis()->GetXmax()); } + if (fData[0].diffFourierPhaseOptReal != 0) { + name = TString(fData[0].diffFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.diffFourierPhaseOptReal = new TH1F(name, name, fData[0].diffFourierPhaseOptReal->GetNbinsX(), + fData[0].diffFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].diffFourierPhaseOptReal->GetXaxis()->GetXmax()); + } // calculate all the average data sets double dval; @@ -3779,6 +3917,20 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierPhase->SetMarkerSize(fData[0].dataFourierPhase->GetMarkerSize()); fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); } + if (fDataAvg.dataFourierPhaseOptReal != 0) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierPhaseOptReal->SetMarkerColor(fData[0].dataFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.dataFourierPhaseOptReal->SetLineColor(fData[0].dataFourierPhaseOptReal->GetLineColor()); + fDataAvg.dataFourierPhaseOptReal->SetMarkerSize(fData[0].dataFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.dataFourierPhaseOptReal->SetMarkerStyle(fData[0].dataFourierPhaseOptReal->GetMarkerStyle()); + } if (fDataAvg.theory != 0) { for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; @@ -3845,6 +3997,20 @@ void PMusrCanvas::HandleAverage() fDataAvg.theoryFourierPhase->SetMarkerSize(fData[0].theoryFourierPhase->GetMarkerSize()); fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle()); } + if (fDataAvg.theoryFourierPhaseOptReal != 0) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierPhaseOptReal->SetMarkerColor(fData[0].theoryFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.theoryFourierPhaseOptReal->SetLineColor(fData[0].theoryFourierPhaseOptReal->GetLineColor()); + fDataAvg.theoryFourierPhaseOptReal->SetMarkerSize(fData[0].theoryFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.theoryFourierPhaseOptReal->SetMarkerStyle(fData[0].theoryFourierPhaseOptReal->GetMarkerStyle()); + } if (fDataAvg.diff != 0) { for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; @@ -3910,83 +4076,25 @@ void PMusrCanvas::HandleAverage() fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size()); } // set marker color, line color, maker size, marker type - fDataAvg.diffFourierPhase->SetMarkerColor(fData[0].dataFourierRe->GetMarkerColor()); - fDataAvg.diffFourierPhase->SetLineColor(fData[0].dataFourierRe->GetLineColor()); - fDataAvg.diffFourierPhase->SetMarkerSize(fData[0].dataFourierRe->GetMarkerSize()); - fDataAvg.diffFourierPhase->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle()); + fDataAvg.diffFourierPhase->SetMarkerColor(fData[0].dataFourierPhase->GetMarkerColor()); + fDataAvg.diffFourierPhase->SetLineColor(fData[0].dataFourierPhase->GetLineColor()); + fDataAvg.diffFourierPhase->SetMarkerSize(fData[0].dataFourierPhase->GetMarkerSize()); + fDataAvg.diffFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); } -} - -//-------------------------------------------------------------------------- -// FindOptimalFourierPhase (private) -//-------------------------------------------------------------------------- -/** - *

The idea to estimate the optimal phase is that the imaginary part of the fourier should be - * an antisymmetric function around the resonance, hence the asymmetry defined as asymmetry = max+min, - * where max/min is the maximum and minimum of the imaginary part, should be a minimum for the correct phase. - */ -double PMusrCanvas::FindOptimalFourierPhase() -{ - // check that Fourier is really present - if ((fData[0].dataFourierRe == 0) || (fData[0].dataFourierIm == 0)) - return 0.0; - - Double_t minPhase, x, valIm, val_xMin = 0.0, val_xMax = 0.0; - Double_t minIm = 0.0, maxIm = 0.0, asymmetry; - // get min/max of the imaginary part for phase = 0.0 as a starting point - minPhase = 0.0; - Bool_t first = true; - for (Int_t i=0; iGetNbinsX(); i++) { - x = fData[0].dataFourierIm->GetBinCenter(i); - if ((x > fFourier.fRangeForPhaseCorrection[0]) && (x < fFourier.fRangeForPhaseCorrection[1])) { - valIm = fData[0].dataFourierIm->GetBinContent(i); - if (first) { - minIm = valIm; - maxIm = valIm; - val_xMin = valIm; - first = false; - } else { - if (valIm < minIm) - minIm = valIm; - if (valIm > maxIm) - maxIm = valIm; - val_xMax = valIm; + if (fDataAvg.diffFourierPhaseOptReal != 0) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); } + fDataAvg.diffFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierPhaseOptReal->SetMarkerColor(fData[0].dataFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.diffFourierPhaseOptReal->SetLineColor(fData[0].dataFourierPhaseOptReal->GetLineColor()); + fDataAvg.diffFourierPhaseOptReal->SetMarkerSize(fData[0].dataFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.diffFourierPhaseOptReal->SetMarkerStyle(fData[0].dataFourierPhaseOptReal->GetMarkerStyle()); } - asymmetry = (maxIm+minIm)*(val_xMin-val_xMax); - - // go through all phases an check if there is a larger max-min value of the imaginary part - double cp, sp; - for (double phase=0.1; phase < 180.0; phase += 0.1) { - cp = TMath::Cos(phase / 180.0 * TMath::Pi()); - sp = TMath::Sin(phase / 180.0 * TMath::Pi()); - first = true; - for (Int_t i=0; iGetNbinsX(); i++) { - x = fData[0].dataFourierIm->GetBinCenter(i); - if ((x > fFourier.fRangeForPhaseCorrection[0]) && (x < fFourier.fRangeForPhaseCorrection[1])) { - valIm = -sp * fData[0].dataFourierRe->GetBinContent(i) + cp * fData[0].dataFourierIm->GetBinContent(i); - if (first) { - minIm = valIm; - maxIm = valIm; - val_xMin = valIm; - first = false; - } else { - if (valIm < minIm) - minIm = valIm; - if (valIm > maxIm) - maxIm = valIm; - val_xMax = valIm; - } - } - } - if (fabs(asymmetry) > fabs((maxIm+minIm)*(val_xMin-val_xMax))) { - minPhase = phase; - asymmetry = (maxIm+minIm)*(val_xMin-val_xMax); - } - } - - return minPhase; } //-------------------------------------------------------------------------- @@ -4030,6 +4138,10 @@ void PMusrCanvas::CleanupFourier() delete fData[i].dataFourierPhase; fData[i].dataFourierPhase = 0; } + if (fData[i].dataFourierPhaseOptReal != 0) { + delete fData[i].dataFourierPhaseOptReal; + fData[i].dataFourierPhaseOptReal = 0; + } if (fData[i].theoryFourierRe != 0) { delete fData[i].theoryFourierRe; fData[i].theoryFourierRe = 0; @@ -4046,6 +4158,10 @@ void PMusrCanvas::CleanupFourier() delete fData[i].theoryFourierPhase; fData[i].theoryFourierPhase = 0; } + if (fData[i].theoryFourierPhaseOptReal != 0) { + delete fData[i].theoryFourierPhaseOptReal; + fData[i].theoryFourierPhaseOptReal = 0; + } } } @@ -4105,6 +4221,10 @@ void PMusrCanvas::CleanupAverage() delete fDataAvg.dataFourierPhase; fDataAvg.dataFourierPhase = 0; } + if (fDataAvg.dataFourierPhaseOptReal != 0) { + delete fDataAvg.dataFourierPhaseOptReal; + fDataAvg.dataFourierPhaseOptReal = 0; + } if (fDataAvg.theory != 0) { delete fDataAvg.theory; fDataAvg.theory = 0; @@ -4125,6 +4245,10 @@ void PMusrCanvas::CleanupAverage() delete fDataAvg.theoryFourierPhase; fDataAvg.theoryFourierPhase = 0; } + if (fDataAvg.theoryFourierPhaseOptReal != 0) { + delete fDataAvg.theoryFourierPhaseOptReal; + fDataAvg.theoryFourierPhaseOptReal = 0; + } if (fDataAvg.diff != 0) { delete fDataAvg.diff; fDataAvg.diff = 0; @@ -4145,6 +4269,10 @@ void PMusrCanvas::CleanupAverage() delete fDataAvg.diffFourierPhase; fDataAvg.diffFourierPhase = 0; } + if (fDataAvg.diffFourierPhaseOptReal != 0) { + delete fDataAvg.diffFourierPhaseOptReal; + fDataAvg.diffFourierPhaseOptReal = 0; + } } //-------------------------------------------------------------------------- @@ -5380,6 +5508,83 @@ void PMusrCanvas::PlotFourier(Bool_t unzoom) } break; + case PV_FOURIER_PHASE_OPT_REAL: + // set x-range + if ((fFourier.fPlotRange[0] != -1) && (fFourier.fPlotRange[1] != -1)) { + xmin = fFourier.fPlotRange[0]; + xmax = fFourier.fPlotRange[1]; + } else { + xmin = fData[0].dataFourierPhaseOptReal->GetBinLowEdge(1); + xmax = fData[0].dataFourierPhaseOptReal->GetBinLowEdge(fData[0].dataFourierPhaseOptReal->GetNbinsX())+fData[0].dataFourierPhaseOptReal->GetBinWidth(1); + } + + // set y-range + // first find minimum/maximum of all histos + ymin = GetMinimum(fData[0].dataFourierPhaseOptReal); + ymax = GetMaximum(fData[0].dataFourierPhaseOptReal); + binContent = GetMinimum(fData[0].theoryFourierPhaseOptReal); + if (binContent < ymin) + ymin = binContent; + binContent = GetMaximum(fData[0].theoryFourierPhaseOptReal); + if (binContent > ymax) + ymax = binContent; + for (UInt_t i=1; i ymax) + ymax = binContent; + binContent = GetMinimum(fData[i].theoryFourierPhaseOptReal); + if (binContent < ymin) + ymin = binContent; + binContent = GetMaximum(fData[i].theoryFourierPhaseOptReal); + if (binContent > ymax) + ymax = binContent; + } + + // delete old fHistoFrame if present + if (fHistoFrame) { + delete fHistoFrame; + fHistoFrame = 0; + } + + fHistoFrame = fDataTheoryPad->DrawFrame(xmin, 1.05*ymin, xmax, 1.05*ymax); + + // find the maximal number of points present in the histograms and increase the default number of points of fHistoFrame (1000) to the needed one + noOfPoints = 1000; + for (UInt_t i=0; iGetNbinsX() > (Int_t)noOfPoints) + noOfPoints = fData[i].dataFourierPhaseOptReal->GetNbinsX(); + } + noOfPoints *= 2; // make sure that there are enough points + fHistoFrame->SetBins(noOfPoints, xmin, xmax); + + for (UInt_t i=0; iGetXaxis()->SetRangeUser(xmin, xmax); + fData[i].dataFourierPhaseOptReal->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + fData[i].theoryFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax); + fData[i].theoryFourierPhaseOptReal->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + } + + // set x-axis title + fHistoFrame->GetXaxis()->SetTitle(xAxisTitle.Data()); + + // set y-axis title + fHistoFrame->GetYaxis()->SetTitleOffset(1.3); + fHistoFrame->GetYaxis()->SetTitle("Phase Opt. Real Fourier"); + + // plot data + for (UInt_t i=0; iDraw("psame"); + } + + // plot theories + for (UInt_t i=0; iDraw("same"); + } + + break; default: break; } @@ -5848,6 +6053,9 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom) case PV_FOURIER_PHASE: yAxisTitle = ""; break; + case PV_FOURIER_PHASE_OPT_REAL: + yAxisTitle = ""; + break; default: yAxisTitle = "??"; break; @@ -5967,6 +6175,14 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom) fDataAvg.diffFourierPhase->Draw("psame"); } break; + case PV_FOURIER_PHASE_OPT_REAL: + if (!fDifferenceView) { // averaged Fourier Phase Opt Real view + fDataAvg.dataFourierPhaseOptReal->Draw("psame"); + fDataAvg.theoryFourierPhaseOptReal->Draw("same"); + } else { // averaged diff Fourier Phase view + fDataAvg.diffFourierPhaseOptReal->Draw("psame"); + } + break; default: break; } diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 42e73dcc..cdc6791a 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -55,12 +55,13 @@ #define XTHEO 0.75 // Current Plot Views -#define PV_DATA 1 -#define PV_FOURIER_REAL 2 -#define PV_FOURIER_IMAG 3 -#define PV_FOURIER_REAL_AND_IMAG 4 -#define PV_FOURIER_PWR 5 -#define PV_FOURIER_PHASE 6 +#define PV_DATA 1 +#define PV_FOURIER_REAL 2 +#define PV_FOURIER_IMAG 3 +#define PV_FOURIER_REAL_AND_IMAG 4 +#define PV_FOURIER_PWR 5 +#define PV_FOURIER_PHASE 6 +#define PV_FOURIER_PHASE_OPT_REAL 7 // Canvas menu id's #define P_MENU_ID_DATA 10001 @@ -71,13 +72,14 @@ #define P_MENU_PLOT_OFFSET 1000 -#define P_MENU_ID_FOURIER_REAL 100 -#define P_MENU_ID_FOURIER_IMAG 101 -#define P_MENU_ID_FOURIER_REAL_AND_IMAG 102 -#define P_MENU_ID_FOURIER_PWR 103 -#define P_MENU_ID_FOURIER_PHASE 104 -#define P_MENU_ID_FOURIER_PHASE_PLUS 105 -#define P_MENU_ID_FOURIER_PHASE_MINUS 106 +#define P_MENU_ID_FOURIER_REAL 100 +#define P_MENU_ID_FOURIER_IMAG 101 +#define P_MENU_ID_FOURIER_REAL_AND_IMAG 102 +#define P_MENU_ID_FOURIER_PWR 103 +#define P_MENU_ID_FOURIER_PHASE 104 +#define P_MENU_ID_FOURIER_PHASE_OPT_REAL 105 +#define P_MENU_ID_FOURIER_PHASE_PLUS 106 +#define P_MENU_ID_FOURIER_PHASE_MINUS 107 //------------------------------------------------------------------------ /** @@ -122,16 +124,19 @@ typedef struct { TH1F *dataFourierIm; ///< imaginary part of the Fourier transform of the data histogram TH1F *dataFourierPwr; ///< power spectrum of the Fourier transform of the data histogram TH1F *dataFourierPhase; ///< phase spectrum of the Fourier transform of the data histogram + TH1F *dataFourierPhaseOptReal; ///< phase optimized real part spectrum Fourier transform of the data histogram TH1F *theory; ///< theory histogram belonging to the data histogram TH1F *theoryFourierRe; ///< real part of the Fourier transform of the theory histogram TH1F *theoryFourierIm; ///< imaginary part of the Fourier transform of the theory histogram TH1F *theoryFourierPwr; ///< power spectrum of the Fourier transform of the theory histogram TH1F *theoryFourierPhase; ///< phase spectrum of the Fourier transform of the theory histogram + TH1F *theoryFourierPhaseOptReal; ///< phase optimized real part spectrum Fourier transform of the theory histogram TH1F *diff; ///< difference histogram, i.e. data-theory TH1F *diffFourierRe; ///< real part of the Fourier transform of the diff histogram TH1F *diffFourierIm; ///< imaginary part of the Fourier transform of the diff histogram TH1F *diffFourierPwr; ///< power spectrum of the Fourier transform of the diff histogram TH1F *diffFourierPhase; ///< phase spectrum of the Fourier transform of the diff histogram + TH1F *diffFourierPhaseOptReal; ///< phase optimized real part spectrum Fourier transform of the diff histogram PMusrCanvasPlotRange *dataRange; ///< keep the msr-file plot data range UInt_t diffFourierTag; ///< 0=not relevant, 1=d-f (Fourier of difference time spectra), 2=f-d (difference of Fourier spectra) } PMusrCanvasDataSet; @@ -302,7 +307,6 @@ class PMusrCanvas : public TObject, public TQObject virtual void HandleDifferenceFourier(); virtual void HandleFourierDifference(); virtual void HandleAverage(); - virtual Double_t FindOptimalFourierPhase(); virtual void CleanupDifference(); virtual void CleanupFourier(); virtual void CleanupFourierDifference(); From 9c1144fa7bf66abe8d3dae8bd522caf0f5de205f Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Wed, 14 Dec 2016 13:57:11 +0100 Subject: [PATCH 06/13] Deal with old pta runs. --- src/external/MuSRFitGUI/MSR.pm | 38 ++++++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/src/external/MuSRFitGUI/MSR.pm b/src/external/MuSRFitGUI/MSR.pm index a6585671..713abcaf 100755 --- a/src/external/MuSRFitGUI/MSR.pm +++ b/src/external/MuSRFitGUI/MSR.pm @@ -446,7 +446,7 @@ $logxy"; if ($All{"FPLOT"} eq $EMPTY) {$All{"FPLOT"}="POWER";} if ($All{"FPHASE"} eq $EMPTY) {$All{"FPHASE"}="8.5";} my $FrqRange = "#range ".$All{"FRQMIN"}." ".$All{"FRQMAX"}; - if ($All{"FRQMAX"} ne $EMPTY && $All{"FRQMIN"} ne $EMPTY) { + if ($All{"FRQMAX"} ne $EMPTY && $All{"FRQMIN"} ne $EMPTY && $All{"FRQMAX"} ne $All{"FRQMIN"}) { $FrqRange = "range ".$All{"FRQMIN"}." ".$All{"FRQMAX"}; } @@ -1780,20 +1780,40 @@ sub RUNFileNameAuto { $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; } elsif ( $BeamLine eq "GPS" ) { - $RUN_File_Name = "deltat_tdc_gps_" . $RUNtmp; - $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + if ($YEAR >= 2008) { + $RUN_File_Name = "deltat_tdc_gps_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + } else { + $RUN_File_Name = "deltat_pta_gps_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/$RUN_File_Name"; + } } elsif ( $BeamLine eq "LTF" ) { - $RUN_File_Name = "deltat_tdc_ltf_" . $RUNtmp; - $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + if ($YEAR >= 2009) { + $RUN_File_Name = "deltat_tdc_ltf_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + } else { + $RUN_File_Name = "deltat_pta_ltf_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/$RUN_File_Name"; + } } elsif ( $BeamLine eq "Dolly" ) { - $RUN_File_Name = "deltat_tdc_dolly_" . $RUNtmp; - $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + if ($YEAR >= 2008) { + $RUN_File_Name = "deltat_tdc_dolly_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + } else { + $RUN_File_Name = "deltat_pta_dolly_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/$RUN_File_Name"; + } } elsif ( $BeamLine eq "GPD" ) { - $RUN_File_Name = "deltat_tdc_gpd_" . $RUNtmp; - $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + if ($YEAR >= 2008) { + $RUN_File_Name = "deltat_tdc_gpd_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/tdc/$RUN_File_Name"; + } else { + $RUN_File_Name = "deltat_pta_gpd_" . $RUNtmp; + $RUNFILE = "$DATADIR/d$YEAR/$RUN_File_Name"; + } } elsif ( $BeamLine eq "HAL" ) { $RUNtmp=sprintf("%05d",$RUNtmp); From 65d40cfe978deb3a3ebd9dbe183b7d928b5ded21 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 18 Dec 2016 10:31:16 +0100 Subject: [PATCH 07/13] improved user interface. (i) it is searching now not only in the current directory or file-path-name, but also in the default search paths for data files as defined in the musrfit_startup.xml. (ii) it is now also possible to use run numbers rather than full file names. --- src/dump_header.cpp | 417 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 392 insertions(+), 25 deletions(-) diff --git a/src/dump_header.cpp b/src/dump_header.cpp index 20f816ea..f6affbac 100644 --- a/src/dump_header.cpp +++ b/src/dump_header.cpp @@ -31,12 +31,16 @@ #include "config.h" #endif +#include +#include +#include #include #include #include #include #include +#include #include #include @@ -46,6 +50,7 @@ using namespace std; #include #include "git-revision.h" +#include "PStartupHandler.h" #include "TMusrRunHeader.h" #include "TLemRunHeader.h" #include "MuSR_td_PSI_bin.h" @@ -64,16 +69,29 @@ using namespace std; */ void dump_header_syntax() { - cout << endl << "usage: dump_header [--file_format ] | --help | --version"; - cout << endl << " Dumps the header information of a given muSR data file onto the standard output."; - cout << endl << " If no is given, it will try to fiddle out what it might be."; + cout << endl << "usage: dump_header [-rn | -fn ] [-ff, --fileFormat ]"; + cout << endl << " [-y, --year ] [-s, --summary] [--psi-bulk ] |"; + cout << endl << " --help | --version"; cout << endl; - cout << endl << " : muSR data file name."; - cout << endl << " --file_format : where can be:"; + cout << endl << " Dumps the header information of a given muSR data file onto the standard output."; + cout << endl << " If no info is povided, it will try to guess what it might be."; + cout << endl << " For guessing of the file format is not possible. The default assumption here is 'MusrRoot'."; + cout << endl; + cout << endl << " -rn, --runNo : run number of the header to be dumped."; + cout << endl << " -fn, --fileName : muSR data file name."; + cout << endl << " -ff, --fileFormat : where can be:"; cout << endl << " MusrRoot, NeXus, ROOT (old LEM), PSI-BIN, PSI-MDU, MUD, WKM"; cout << endl << " NeXus is only supported if enabled."; - cout << endl << " --help, -h : will show this help"; - cout << endl << " --version, -v : will show the current version."; + cout << endl << " -y, --year : has to be 4 digit, e.g. 2005, if provided it is used to"; + cout << endl << " generate the file name for the given , otherwise the current"; + cout << endl << " year is used. If a file name is given, this option has no effect."; + cout << endl << " -s, --summary : this option is used for LE-uSR data sets only. It will, additionally"; + cout << endl << " to the header information, print the summary file content."; + cout << endl << " --psi-bulk : where consists of two items: (i) pta or tdc, "; + cout << endl << " (ii) gps | ltf | dolly | gpd | hifi. This is needed in combination with"; + cout << endl << " the file formats PSI-BIN and PSI-MDU."; + cout << endl << " -h, --help : will show this help"; + cout << endl << " -v, --version : will show the current version."; cout << endl << endl; } @@ -81,7 +99,7 @@ void dump_header_syntax() /** * */ -int dump_header_root(const string fileName, const string fileFormat) +int dump_header_root(const string fileName, const string fileFormat, const bool summary) { TFile f(fileName.c_str()); if (f.IsZombie()) { @@ -178,6 +196,27 @@ int dump_header_root(const string fileName, const string fileFormat) delete header; } + // summary as well? + if (summary && (fileType == DH_MUSR_ROOT)) { + TObjArray *runSum=0; + runSum = (TObjArray*)folder->FindObject("RunSummary"); + if (!runSum) { // something is wrong!! + cerr << endl << "**ERROR** Couldn't obtain RunSummary " << fileName << endl; + f.Close(); + return 1; + } + cout << "++++++++++++++++++++" << endl; + cout << " Run Summary" << endl; + cout << "++++++++++++++++++++" << endl; + TObjString *tstr; + TString str; + for (Int_t i=0; iGetEntries(); i++) { + tstr = (TObjString*)runSum->At(i); + str = tstr->String(); + cout << str; + } + } + f.Close(); return 0; @@ -562,7 +601,145 @@ int dump_header_wkm(const string fileName, const string fileFormat) //------------------------------------------------------------------------ /** - * + * @brief is_number + * @param s + * @return + */ +bool dump_is_number(const char *s) +{ + int i=0; + + if (s == 0) // make sure it is not a null pointer + return false; + + while (isdigit(s[i])) + i++; + + if (s[i] == '\0') + return true; + else + return false; +} + +//------------------------------------------------------------------------ +/** + * @brief dump_current_year + * @return + */ +int dump_current_year() +{ + time_t rawtime; + struct tm * timeinfo; + char buffer[32]; + + time (&rawtime); + timeinfo = localtime(&rawtime); + strftime(buffer, 32, "%Y", timeinfo); + + return atoi(buffer); +} + +//------------------------------------------------------------------------ +/** + * @brief dump_create_fln + * @param runNo + * @param year + * @param fileFormat + * @return + */ +string dump_create_fln(string runNo, string year, string fileFormat, bool pta, string instrument) +{ + string result = "??"; + int yearShort=0; + int iRunNo=0; + + if (fileFormat.empty()) + fileFormat = "MusrRoot"; + + // make sure that a 'legal' file format has been found + if (!boost::iequals(fileFormat, "MusrRoot") && + !boost::iequals(fileFormat, "NeXus") && + !boost::iequals(fileFormat, "ROOT") && + !boost::iequals(fileFormat, "PSI-BIN") && + !boost::iequals(fileFormat, "PSI-MDU") && + !boost::iequals(fileFormat, "MDU") && + !boost::iequals(fileFormat, "WKM")) { + return result; + } + + // if year is an empty string get the current year + int yy=-1; + stringstream ss; + if (year.empty()) { + yy = dump_current_year(); + ss << yy; + year = ss.str(); + } + yy = atoi(year.c_str()); + if (yy > 2000) + yearShort = yy - 2000; + else + yearShort = yy - 1900; + + iRunNo = atoi(runNo.c_str()); + + char fln[64]; + char ptatdc[8]; + if (pta) + strcpy(ptatdc, "pta"); + else + strcpy(ptatdc, "tdc"); + if (boost::iequals(fileFormat, "MusrRoot") || boost::iequals(fileFormat, "ROOT")) { + snprintf(fln, sizeof(fln), "lem%02d_his_%04d.root", yearShort, iRunNo); + } else if (boost::iequals(fileFormat, "NeXus")) { + snprintf(fln, sizeof(fln), "%s.nxs", runNo.c_str()); + } else if (boost::iequals(fileFormat, "PSI-BIN")) { + snprintf(fln, sizeof(fln), "deltat_%s_%s_%04d.bin", ptatdc, instrument.c_str(), iRunNo); + } else if (boost::iequals(fileFormat, "PSI-MDU")) { + snprintf(fln, sizeof(fln), "%s_%s_%s_%05d.mdu", ptatdc, instrument.c_str(), year.c_str(), iRunNo); + } else if (boost::iequals(fileFormat, "MUD")) { + snprintf(fln, sizeof(fln), "%06d.msr", iRunNo); + } else if (boost::iequals(fileFormat, "WKM")) { + + } + result = fln; + + return result; +} + +//------------------------------------------------------------------------ +/** + * @brief dump_file_exists + * @param pathName + * @return + */ +bool dump_file_exists(const string pathName) +{ + bool exists = true; + + int res = access(pathName.c_str(), R_OK); + if (res < 0) { + if (errno == ENOENT) { + // file does not exist + exists = false; + } else if (errno == EACCES) { + // file exists but is not readable + exists = false; + } else { + // FAIL + exists = false; + } + } + + return exists; +} + +//------------------------------------------------------------------------ +/** + * @brief main + * @param argc + * @param argv + * @return */ int main(int argc, char *argv[]) { @@ -571,9 +748,14 @@ int main(int argc, char *argv[]) return 0; } + string runNo(""); string fileName(""); string fileFormat(""); - int count=0; + string year(""); + bool pta(false); + string instrument(""); + bool summary(false); + for (int i=1; i= argc) { - cout << endl << "**ERROR** found '--file_format' without !" << endl; + cerr << endl << "**ERROR** found -rn, --runNo without !" << endl; + dump_header_syntax(); + return 1; + } + // make sure there is only one 'legal' run number + int count = 1; + while (dump_is_number(argv[i+count]) && (i+count < argc)) + count++; + // make sure there is one and only one run number given + if (count == 1) { + cerr << endl << "**ERROR** found -rn, --runNo without , or the provided ('" << argv[i+1] << "') is not a number!" << endl; + dump_header_syntax(); + return 1; + } + if (count > 2) { + cerr << endl << "**ERROR** found -rn, --runNo with more than one ! This is not yet supported." << endl; + dump_header_syntax(); + return 1; + } + runNo = argv[i+1]; + i++; + } else if (!strcmp(argv[i], "-fn") || !strcmp(argv[i], "--fileName")) { + if (i+1 >= argc) { + cerr << endl << "**ERROR** found -fn, --fileName without !" << endl; + dump_header_syntax(); + return 1; + } + fileName = argv[i+1]; + i++; + } else if (!strcmp(argv[i], "--fileFormat") || !strcmp(argv[i], "-ff")) { + if (i+1 >= argc) { + cerr << endl << "**ERROR** found -ff, --fileFormat without !" << endl; dump_header_syntax(); return 1; } @@ -595,27 +808,169 @@ int main(int argc, char *argv[]) if (!boost::iequals(ff, "MusrRoot") && !boost::iequals(ff, "NeXus") && !boost::iequals(ff, "ROOT") && !boost::iequals(ff, "PSI-BIN") && !boost::iequals(ff, "PSI-MDU") && !boost::iequals(ff, "MUD") && !boost::iequals(ff, "WKM")) { // none of the listed found - cout << endl << "**ERROR** found unsupported muSR file data format: " << argv[i+1] << endl; + cerr << endl << "**ERROR** found unsupported muSR file data format: " << argv[i+1] << endl; dump_header_syntax(); return 1; } fileFormat = argv[i+1]; i++; + } else if (!strcmp(argv[i], "-y") || !strcmp(argv[i], "--year")) { + if (i+1 >= argc) { + cerr << endl << "**ERROR** found -y, --year without !" << endl; + dump_header_syntax(); + return 1; + } + if (!dump_is_number(argv[i+1])) { + cerr << endl << "**ERROR** found -y, --year with sensless '" << argv[i+1] << "'!" << endl; + dump_header_syntax(); + return 1; + } + int yy = strtod(argv[i+1], (char**)0); + if ((yy < 1950) || (yy > dump_current_year())) { + cerr << endl << "**ERROR** found -y, --year with '" << yy << "'!"; + cerr << endl << " Well, cannot handle files in the pre-muSR time nor in the future." << endl; + dump_header_syntax(); + return 1; + } + year = argv[i+1]; + i++; + } else if (!strcmp(argv[i], "-s") || !strcmp(argv[i], "--summary")) { + summary = true; + } else if (!strcmp(argv[i], "--psi-bulk")) { + if (i+2 >= argc) { + cerr << endl << "**ERROR** found --psi-bulk with insufficient input!" << endl; + dump_header_syntax(); + return 1; + } + if (!strcmp(argv[i+1], "pta")) + pta = true; + else if (!strcmp(argv[i+1], "tdc")) + pta = false; + else { + cerr << endl << "**ERROR** found --psi-bulk with 1st argument '" << argv[i+1] << "'! Allowed is 'pta' or 'tdc'." << endl; + dump_header_syntax(); + return 1; + } + if (strcmp(argv[i+2], "gps") && strcmp(argv[i+2], "ltf") && strcmp(argv[i+2], "dolly") && + strcmp(argv[i+2], "gpd") && strcmp(argv[i+2], "hifi")) { + cerr << endl << "**ERROR** found --psi-bulk with 2nd argument '" << argv[i+1] << "'! This is an unkown instrument." << endl; + dump_header_syntax(); + return 1; + } + instrument = argv[i+2]; + i += 2; } else { - count++; - fileName = argv[i]; + cerr << endl << "**ERROR** found unkown option '" << argv[i] << "'." << endl; + dump_header_syntax(); + return 1; } } - // check if more then one file name was given - if (count != 1) { - cout << endl << "**ERROR** (only) a single file name is needed!" << endl; - dump_header_syntax(); + // if year is not provided, take the current one + if (year.empty()) { + stringstream ss; + ss << dump_current_year(); + year = ss.str(); + } + + if ((runNo.length() != 0) && (fileName.length() != 0)) { + cerr << endl << "**ERROR** currently only either runNo or fileName can be handled, not both simultanously." << endl; + return 1; + } + + // invoke the startup handler in order to get the default search paths to the data files + // read startup file + char startup_path_name[128]; + TSAXParser *saxParser = new TSAXParser(); + PStartupHandler *startupHandler = new PStartupHandler(); + if (!startupHandler->StartupFileFound()) { + cerr << endl << ">> musrfit **WARNING** couldn't find " << startupHandler->GetStartupFilePath().Data(); + cerr << endl; + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } + } else { + strcpy(startup_path_name, startupHandler->GetStartupFilePath().Data()); + saxParser->ConnectToHandler("PStartupHandler", startupHandler); + //status = saxParser->ParseFile(startup_path_name); + // parsing the file as above seems to lead to problems in certain environments; + // use the parseXmlFile function instead (see PStartupHandler.cpp for the definition) + int status = parseXmlFile(saxParser, startup_path_name); + // check for parse errors + if (status) { // error + cerr << endl << ">> musrfit **WARNING** Reading/parsing musrfit_startup.xml failed."; + cerr << endl; + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } + } + } + + // runNo given, hence try to create the necessary file name based on the provided information + if (runNo != "") { + string str = dump_create_fln(runNo, year, fileFormat, pta, instrument); + if (str == "??") { + cerr << endl << "**ERROR** couldn't get a proper file name." << endl; + return 1; + } + fileName = str; + } + + bool found_fln = false; + // 1st check if the file name is found in the current directory + string pathFln(""); + pathFln = "./" + fileName; + if (dump_file_exists(pathFln)) + found_fln = true; + + // 2nd check if file name is found in any default search paths if not already found in the current directory + if (!found_fln) { + PStringVector pathList = startupHandler->GetDataPathList(); + for (unsigned int i=0; i Date: Sun, 18 Dec 2016 10:36:02 +0100 Subject: [PATCH 08/13] for maxLH it is now possible to get the per-run-block maxLH, and the expected maxLH --- src/classes/PFitter.cpp | 104 +++++----- src/classes/PFitterFcn.cpp | 11 +- src/classes/PMsrHandler.cpp | 305 ++++++++++++++++------------- src/classes/PRunListCollection.cpp | 93 ++++++++- src/classes/PRunSingleHisto.cpp | 100 +++++++++- src/include/PMsrHandler.h | 2 +- src/include/PMusrCanvas.h | 6 +- src/include/PRunListCollection.h | 3 + src/include/PRunSingleHisto.h | 1 + 9 files changed, 433 insertions(+), 192 deletions(-) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index e366bda1..c921c250 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -1564,64 +1564,66 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) } // handle expected chisq if applicable - if (fUseChi2) { - fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back + fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back - // calculate expected chisq - std::vector param; - Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerHisto; - std::vector ndfPerHisto; + // calculate expected chisq + std::vector param; + Double_t totalExpectedChisq = 0.0; + std::vector expectedChisqPerHisto; + std::vector ndfPerHisto; - for (UInt_t i=0; iCalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + // CalcExpectedChiSquare handles both, chisq and mlh + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); - // calculate chisq per run - std::vector chisqPerHisto; - for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + // calculate chisq per run + std::vector chisqPerHisto; + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + if (fUseChi2) chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i)); - } - - if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only - // get the ndf's of the histos - UInt_t ndf_histo; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - ndfPerHisto.push_back(ndf_histo); - } - - // feed the msr-file handler - PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); - if (statistics) { - statistics->fMinPerHisto = chisqPerHisto; - statistics->fMinExpected = totalExpectedChisq; - statistics->fMinExpectedPerHisto = expectedChisqPerHisto; - statistics->fNdfPerHisto = ndfPerHisto; - } - } else if (chisqPerHisto.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits - UInt_t ndf_histo = 0; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - ndfPerHisto.push_back(ndf_histo); - } - - // feed the msr-file handler - PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); - if (statistics) { - statistics->fMinPerHisto = chisqPerHisto; - statistics->fNdfPerHisto = ndfPerHisto; - } - } - - // clean up - param.clear(); - expectedChisqPerHisto.clear(); - ndfPerHisto.clear(); - chisqPerHisto.clear(); + else + chisqPerHisto.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); } + if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only + // get the ndf's of the histos + UInt_t ndf_histo; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_histo); + } + + // feed the msr-file handler + PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); + if (statistics) { + statistics->fMinPerHisto = chisqPerHisto; + statistics->fMinExpected = totalExpectedChisq; + statistics->fMinExpectedPerHisto = expectedChisqPerHisto; + statistics->fNdfPerHisto = ndfPerHisto; + } + } else if (chisqPerHisto.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_histo); + } + + // feed the msr-file handler + PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); + if (statistics) { + statistics->fMinPerHisto = chisqPerHisto; + statistics->fNdfPerHisto = ndfPerHisto; + } + } + + // clean up + param.clear(); + expectedChisqPerHisto.clear(); + ndfPerHisto.clear(); + chisqPerHisto.clear(); + cout << ">> PFitter::ExecuteSave(): will write minuit2 output file ..." << endl; ofstream fout; diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index c6a8566a..07c4831a 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -110,15 +110,20 @@ void PFitterFcn::CalcExpectedChiSquare(const std::vector &par, Double_ totalExpectedChisq = 0.0; expectedChisqPerRun.clear(); - // only do something for chisq + Double_t value = 0.0; if (fUseChi2) { - Double_t value = 0.0; - // single histo for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i' expectedChisqPerRun.push_back(value); totalExpectedChisq += value; } + } else { // log max. likelihood + // single histo + for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { + value = fRunListCollection->GetSingleHistoMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i' + expectedChisqPerRun.push_back(value); + totalExpectedChisq += value; + } } } diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index d126ccd9..dedf2d29 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -825,13 +825,7 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout << left << "lifetime"; fout << fRuns[runNo].GetLifetimeParamNo() << endl; } else if (sstr.BeginsWith("lifetimecorrection")) { -/* obsolate, hence do nothing here - if ((fRuns[runNo].IsLifetimeCorrected()) && - ((fRuns[runNo].GetFitType() == MSR_FITTYPE_SINGLE_HISTO) || - (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO))) { - fout << "lifetimecorrection" << endl; - } -*/ + // obsolate, hence do nothing here } else if (sstr.BeginsWith("map")) { fout << "map "; for (UInt_t j=0; jsize(); j++) { @@ -1215,38 +1209,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) if (fStatistic.fValid) { // valid fit result if (fStatistic.fChisq) { str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; + } else { + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); + } + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; + + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + if (fStatistic.fChisq) { + str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + } else { + str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + } + if (fStartupOptions) { + if (fStartupOptions->writeExpectedChisq) + fout << str.Data() << endl; + } if (messages) cout << endl << str.Data() << endl; - // check if expected chisq needs to be written - if (fStatistic.fMinExpected != 0.0) { - str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", - fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - if (messages) - cout << endl << str.Data() << endl; - - for (UInt_t i=0; i 0) { + for (UInt_t i=0; i 0) { + if (fStatistic.fChisq) { str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) = (%d/%lf/%lf)", i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) = (%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); } - } - } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written - for (UInt_t i=0; iwriteExpectedChisq) fout << str.Data() << endl; @@ -1256,11 +1250,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } - } else { // maxLH - str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; iwriteExpectedChisq) + fout << str.Data() << endl; + } + + if (messages) + cout << str.Data() << endl; + } } } else { fout << "*** FIT DID NOT CONVERGE ***" << endl; @@ -1272,38 +1278,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) if (fStatistic.fValid) { // valid fit result if (fStatistic.fChisq) { // chisq str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else { + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); + } + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; - // check if expected chisq needs to be written - if (fStatistic.fMinExpected != 0.0) { + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + if (fStatistic.fChisq) { // chisq str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + } + if (fStartupOptions) { + if (fStartupOptions->writeExpectedChisq) + fout << str.Data() << endl; + } + if (messages) + cout << str.Data() << endl; - for (UInt_t i=0; i 0) { + for (UInt_t i=0; i 0) { + if (fStatistic.fChisq) { // chisq str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) = (%d/%lf/%lf)", i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) = (%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); } - } - } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written - for (UInt_t i=0; iwriteExpectedChisq) fout << str.Data() << endl; @@ -1313,11 +1319,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } - } else { // max. log. liklihood - str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; iwriteExpectedChisq) + fout << str.Data() << endl; + } + + if (messages) + cout << str.Data() << endl; + } } } else { fout << "*** FIT DID NOT CONVERGE ***" << endl; @@ -1328,7 +1346,7 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) if (str.Length() > 0) { sstr = str; sstr.Remove(TString::kLeading, ' '); - if (!sstr.BeginsWith("expected chisq") && !sstr.BeginsWith("run block")) + if (!sstr.BeginsWith("expected chisq") && !sstr.BeginsWith("expected maxLH") && !sstr.BeginsWith("run block")) fout << str.Data() << endl; } else { // only write endl if not eof is reached. This is preventing growing msr-files, i.e. more and more empty lines at the end of the file if (!fin.eof()) @@ -1351,38 +1369,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) if (fStatistic.fValid) { // valid fit result if (fStatistic.fChisq) { str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else { + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); + } + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; - // check if expected chisq needs to be written - if (fStatistic.fMinExpected != 0.0) { + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + if (fStatistic.fChisq) { str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + } + if (fStartupOptions) { + if (fStartupOptions->writeExpectedChisq) + fout << str.Data() << endl; + } + if (messages) + cout << str.Data() << endl; - for (UInt_t i=0; i 0) { + for (UInt_t i=0; i 0) { + if (fStatistic.fChisq) { str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) = (%d/%lf/%lf)", i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) = (%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); } - } - } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written - for (UInt_t i=0; iwriteExpectedChisq) fout << str.Data() << endl; @@ -1392,13 +1410,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } - } else { // maxLH - str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; iwriteExpectedChisq) + fout << str.Data() << endl; + } - // check if per run block maxLH needs to be written + if (messages) + cout << str.Data() << endl; + } } } else { fout << "*** FIT DID NOT CONVERGE ***" << endl; @@ -1416,38 +1444,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) if (fStatistic.fValid) { // valid fit result if (fStatistic.fChisq) { // chisq str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else { + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); + } + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; - // check if expected chisq needs to be written - if (fStatistic.fMinExpected != 0.0) { + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + if (fStatistic.fChisq) { // chisq str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + } + if (fStartupOptions) { + if (fStartupOptions->writeExpectedChisq) + fout << str.Data() << endl; + } + if (messages) + cout << str.Data() << endl; - for (UInt_t i=0; i 0) { + for (UInt_t i=0; i 0) { + if (fStatistic.fChisq) { // chisq str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) =(%d/%lf/%lf)", i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); - if (fStartupOptions) { - if (fStartupOptions->writeExpectedChisq) - fout << str.Data() << endl; - } - - if (messages) - cout << str.Data() << endl; + } else { + str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) =(%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); } - } - } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written - for (UInt_t i=0; iwriteExpectedChisq) fout << str.Data() << endl; @@ -1457,11 +1485,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } - } else { // max. log. liklihood - str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; iwriteExpectedChisq) + fout << str.Data() << endl; + } + + if (messages) + cout << str.Data() << endl; + } } } else { fout << "*** FIT DID NOT CONVERGE (4) ***" << endl; @@ -4675,7 +4715,8 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines) if (tstr.Length() > 0) { if (!tstr.BeginsWith("#") && !tstr.BeginsWith("STATISTIC") && !tstr.BeginsWith("chisq") && !tstr.BeginsWith("maxLH") && !tstr.BeginsWith("*** FIT DID NOT CONVERGE ***") && - !tstr.BeginsWith("expected chisq") && !tstr.BeginsWith("run block")) { + !tstr.BeginsWith("expected chisq") && !tstr.BeginsWith("expected maxLH") && + !tstr.BeginsWith("run block")) { cerr << endl << ">> PMsrHandler::HandleStatisticEntry: **SYNTAX ERROR** in line " << lines[i].fLineNo; cerr << endl << ">> '" << lines[i].fLine.Data() << "'"; cerr << endl << ">> not a valid STATISTIC block line"; @@ -6059,6 +6100,7 @@ Bool_t PMsrHandler::EstimateN0() /** *

returns alpha to estimate N0 */ +/*as Double_t PMsrHandler::GetAlphaEstimateN0() { if (fStartupOptions == 0) @@ -6066,6 +6108,7 @@ Double_t PMsrHandler::GetAlphaEstimateN0() return fStartupOptions->alphaEstimateN0; } +*/ //-------------------------------------------------------------------------- // NeededPrecision (private) diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index c9c5e56f..3cc36d01 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -355,7 +355,7 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vectorGetMsrRunList()->at(idx).GetFitType(); - if (type == -1) { // i.e. not forun in the RUN block, try the GLOBAL block + if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->GetFitType(); } @@ -583,6 +583,97 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vectorCalculates expected mlh of the single histogram with run block index idx of a msr-file. + * + * return: + * - expected mlh of for a single histogram + * + * \param par fit parameter vector + * \param idx run block index + */ +Double_t PRunListCollection::GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const +{ + Double_t expected_mlh = 0.0; + + if (idx > fMsrInfo->GetMsrRunList()->size()) { + cerr << ">> PRunListCollection::GetSingleHistoMaximumLikelihoodExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl; + return expected_mlh; + } + + Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); + if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block + type = fMsrInfo->GetMsrGlobal()->GetFitType(); + } + + // count how many entries of this fit-type are present up to idx + UInt_t subIdx = 0; + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } + + // return the mlh of the single run + switch (type) { + case PRUN_SINGLE_HISTO: + expected_mlh = fRunSingleHistoList[subIdx]->CalcMaxLikelihoodExpected(par); + break; + default: + break; + } + + return expected_mlh; +} + +//-------------------------------------------------------------------------- +// GetSingleRunMaximumLikelihood (public) +//-------------------------------------------------------------------------- +/** + *

Calculates mlh of a single run-block entry of the msr-file. + * + * return: + * - mlh of single run-block entry with index idx + * + * \param par fit parameter vector + * \param idx run block index + */ +Double_t PRunListCollection::GetSingleRunMaximumLikelihood(const std::vector& par, const UInt_t idx) const +{ + Double_t mlh = 0.0; + + if (idx > fMsrInfo->GetMsrRunList()->size()) { + cerr << ">> PRunListCollection::GetSingleRunMaximumLikelihood() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl; + return mlh; + } + + Int_t subIdx = 0; + Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); + if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block + type = fMsrInfo->GetMsrGlobal()->GetFitType(); + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } + } + + // return the mlh of the single run + switch (type) { + case PRUN_SINGLE_HISTO: + mlh = fRunSingleHistoList[subIdx]->CalcMaxLikelihood(par); + break; + default: + break; + } + + return mlh; +} + //-------------------------------------------------------------------------- // GetNoOfBinsFitted (public) //-------------------------------------------------------------------------- diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 8b8b3094..23738b70 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -393,6 +393,103 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) return 2.0*mllh; } +//-------------------------------------------------------------------------- +// CalcMaxLikelihoodExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected log maximum-likelihood. + * + * return: + * - log maximum-likelihood value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector& par) +{ + Double_t mllh = 0.0; // maximum log likelihood assuming poisson distribution for the single bin + + Double_t N0; + + // check if norm is a parameter or a function + if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter + N0 = par[fRunInfo->GetNormParamNo()-1]; + } else { // norm is a function + // get function number + UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + + // get tau + Double_t tau; + if (fRunInfo->GetLifetimeParamNo() != -1) + tau = par[fRunInfo->GetLifetimeParamNo()-1]; + else + tau = PMUON_LIFETIME; + + // get background + Double_t bkg; + if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->GetBkgFix(0); + } + } else { // bkg fitted + bkg = par[fRunInfo->GetBkgFitParamNo()-1]; + } + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + Int_t funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par); + } + + // calculate maximum log likelihood + Double_t theo; + Double_t data; + Double_t time(1.0); + Int_t i; + + // norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns + Double_t normalizer = 1.0; + + if (fScaleN0AndBkg) + normalizer = fPacking * (fTimeResolution * 1.0e3); + + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. + // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once + // for a given set of parameters---which should be done outside of the parallelized loop. + // For all other functions it means a tiny and acceptable overhead. + time = fTheory->Func(time, par, fFuncValues); + + #ifdef HAVE_GOMP + Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i,time,theo,data) schedule(dynamic,chunk) reduction(-:mllh) + #endif + for (i=fStartTimeBin; iFunc(time, par, fFuncValues))+bkg; + theo *= normalizer; + + data = normalizer*fData.GetValue()->at(i); + + if (theo <= 0.0) { + cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << endl; + continue; + } + + if (data > 1.0e-9) { // is this correct?? needs to be checked. See G-test + mllh += data*log(data/theo); + } + } + + return 2.0*mllh; +} + //-------------------------------------------------------------------------- // CalcTheory (public) //-------------------------------------------------------------------------- @@ -1571,7 +1668,6 @@ void PRunSingleHisto::EstimateN0() Double_t tau = PMUON_LIFETIME; UInt_t t0 = (UInt_t)round(fT0s[0]); - Double_t alpha = fMsrInfo->GetAlphaEstimateN0(); Double_t dval = 0.0; Double_t nom = 0.0; Double_t denom = 0.0; @@ -1580,14 +1676,12 @@ void PRunSingleHisto::EstimateN0() // calc nominator for (UInt_t i=t0; i 0) denom += xx*xx/dval; diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index 1dab817b..097ccdc1 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -107,7 +107,7 @@ class PMsrHandler virtual void GetGroupingString(Int_t runNo, TString detector, TString &groupingStr); virtual Bool_t EstimateN0(); - virtual Double_t GetAlphaEstimateN0(); +//as virtual Double_t GetAlphaEstimateN0(); private: Bool_t fFourierOnly; ///< flag indicating if Fourier transform only is wished. If yes, some part of the msr-file blocks are not needed. diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index cdc6791a..d61392c2 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -206,12 +206,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 avg=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 avg=false); virtual ~PMusrCanvas(); virtual Bool_t IsValid() { return fValid; } @@ -236,6 +236,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 fStartWithAvg; ///< flag if true, the averaged data/Fourier will be presented 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 @@ -312,6 +313,7 @@ class PMusrCanvas : public TObject, public TQObject virtual void CleanupFourierDifference(); virtual void CleanupAverage(); + virtual void CalcPhaseOptReFT(); virtual Double_t CalculateDiff(const Double_t x, const Double_t y, TH1F *theo); virtual Double_t CalculateDiff(const Double_t x, const Double_t y, TGraphErrors *theo); virtual Int_t FindBin(const Double_t x, TGraphErrors *graph); diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 00b1112e..1c22ee51 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -76,6 +76,9 @@ class PRunListCollection virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; + virtual Double_t GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleRunMaximumLikelihood(const std::vector& par, const UInt_t idx) const; + virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const; virtual UInt_t GetTotalNoOfBinsFitted() const; diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 15423636..ba8b49d4 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -45,6 +45,7 @@ class PRunSingleHisto : public PRunBase virtual Double_t CalcChiSquare(const std::vector& par); virtual Double_t CalcChiSquareExpected(const std::vector& par); virtual Double_t CalcMaxLikelihood(const std::vector& par); + virtual Double_t CalcMaxLikelihoodExpected(const std::vector& par); virtual void CalcTheory(); virtual UInt_t GetNoOfFitBins(); From 63516fc499d3f666bb17bbcdfa8d48186413a988 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 18 Dec 2016 10:37:50 +0100 Subject: [PATCH 09/13] added a first version of a optimized phase correction for the real Fourier transform --- src/classes/PFourier.cpp | 83 ++++++------ src/classes/PFourierCanvas.cpp | 215 +++++++++++++++++++++++-------- src/classes/PMusrCanvas.cpp | 224 +++++++++++++++++++++------------ src/include/PFourier.h | 7 +- src/include/PFourierCanvas.h | 1 + src/musrFT.cpp | 8 +- src/musrview.cpp | 8 +- 7 files changed, 368 insertions(+), 178 deletions(-) diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index b6034b26..8e1ad423 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -94,8 +94,8 @@ void PFTPhaseCorrection::Minimize() // create Minuit2 parameters ROOT::Minuit2::MnUserParameters upar; - upar.Add("c0", fPh_c0, 0.5); - upar.Add("c1", fPh_c1, 0.5); + upar.Add("c0", fPh_c0, 2.0); + upar.Add("c1", fPh_c1, 2.0); // create minimizer ROOT::Minuit2::MnMinimize mn_min(*this, upar); @@ -296,7 +296,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi fValid = true; fIn = 0; fOut = 0; - fPhCorrectedReFT = 0; +//as fPhCorrectedReFT = 0; fApodization = F_APODIZATION_NONE; @@ -384,8 +384,8 @@ PFourier::~PFourier() fftw_free(fIn); if (fOut) fftw_free(fOut); - if (fPhCorrectedReFT) - delete fPhCorrectedReFT; +//as if (fPhCorrectedReFT) +//as delete fPhCorrectedReFT; } //-------------------------------------------------------------------------- @@ -485,73 +485,79 @@ TH1F* PFourier::GetRealFourier(const Double_t scale) } //-------------------------------------------------------------------------- -// GetPhaseOptRealFourier (public) +// GetPhaseOptRealFourier (public, static) //-------------------------------------------------------------------------- /** *

returns the phase corrected real Fourier transform. * * \return the TH1F histo of the phase 'optimzed' real Fourier transform. * + * \param re real part Fourier histogram + * \param im imaginary part Fourier histogram * \param phase return value of the optimal phase dispersion phase[0]+phase[1]*i/N * \param scale normalisation factor * \param min minimal freq / field from which to optimise. Given in the choosen unit. * \param max maximal freq / field up to which to optimise. Given in the choosen unit. */ -TH1F* PFourier::GetPhaseOptRealFourier(vector &phase, const Double_t scale, const Double_t min, const Double_t max) +TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, vector &phase, + const Double_t scale, const Double_t min, const Double_t max) { + if ((re == 0) || (im == 0)) + return 0; + phase.resize(2); // c0, c1 - UInt_t noOfFourierBins = 0; - if (fNoOfBins % 2 == 0) - noOfFourierBins = fNoOfBins/2; - else - noOfFourierBins = (fNoOfBins+1)/2; + TAxis *axis = re->GetXaxis(); - UInt_t minBin = 0; - UInt_t maxBin = noOfFourierBins; + Int_t minBin = 1; + Int_t maxBin = axis->GetNbins(); + Int_t noOfBins = axis->GetNbins(); + Double_t res = axis->GetBinWidth(1); // check if minimum frequency is given. If yes, get the proper minBin if (min != -1.0) { - minBin = (UInt_t)(min/fResolution); + minBin = axis->FindFixBin(min); + if ((minBin == 0) || (minBin > maxBin)) { + minBin = 1; + cerr << "**WARNING** minimum frequency/field out of range. Will adopted it." << endl; + } } // check if maximum frequency is given. If yes, get the proper maxBin if (max != -1.0) { - maxBin = (UInt_t)(max/fResolution); - if (maxBin >= noOfFourierBins) { - maxBin = noOfFourierBins; + maxBin = axis->FindFixBin(max); + if ((maxBin == 0) || (maxBin > axis->GetNbins())) { + maxBin = axis->GetNbins(); cerr << "**WARNING** maximum frequency/field out of range. Will adopted it." << endl; } } // copy the real/imag Fourier from min to max vector realF, imagF; - for (UInt_t i=minBin; i<=maxBin; i++) { - realF.push_back(fOut[i][0]); - imagF.push_back(fOut[i][1]); + for (Int_t i=minBin; i<=maxBin; i++) { + realF.push_back(re->GetBinContent(i)); + imagF.push_back(im->GetBinContent(i)); } // optimize real FT phase - fPhCorrectedReFT = new PFTPhaseCorrection(realF, imagF); - if (fPhCorrectedReFT == 0) { - fValid = false; + PFTPhaseCorrection *phCorrectedReFT = new PFTPhaseCorrection(realF, imagF); + if (phCorrectedReFT == 0) { cerr << endl << "**SEVERE ERROR** couldn't invoke PFTPhaseCorrection object." << endl; return 0; } - fPhCorrectedReFT->Minimize(); - if (!fPhCorrectedReFT->IsValid()) { - fValid = false; - cerr << endl << "**ERROR** could not fina a valid phase correction minimum." << endl; + phCorrectedReFT->Minimize(); + if (!phCorrectedReFT->IsValid()) { + cerr << endl << "**ERROR** could not find a valid phase correction minimum." << endl; return 0; } - phase[0] = fPhCorrectedReFT->GetPhaseCorrectionParam(0); - phase[1] = fPhCorrectedReFT->GetPhaseCorrectionParam(1); + phase[0] = phCorrectedReFT->GetPhaseCorrectionParam(0); + phase[1] = phCorrectedReFT->GetPhaseCorrectionParam(1); // clean up - if (fPhCorrectedReFT) { - delete fPhCorrectedReFT; - fPhCorrectedReFT = 0; + if (phCorrectedReFT) { + delete phCorrectedReFT; + phCorrectedReFT = 0; } realF.clear(); imagF.clear(); @@ -559,21 +565,20 @@ TH1F* PFourier::GetPhaseOptRealFourier(vector &phase, const Double_t s // invoke the real phase optimised histo to be filled. Caller is the owner! Char_t name[256]; Char_t title[256]; - snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", fData->GetName()); - snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", fData->GetTitle()); + snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", re->GetName()); + snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", re->GetTitle()); - TH1F *realPhaseOptFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0); + TH1F *realPhaseOptFourier = new TH1F(name, title, noOfBins, -res/2.0, (Double_t)(noOfBins-1)*res+res/2.0); if (realPhaseOptFourier == 0) { - fValid = false; cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << endl; return 0; } // fill realFourier vector Double_t ph; - for (UInt_t i=0; iSetBinContent(i+1, scale*(fOut[i][0]*cos(ph) - fOut[i][1]*sin(ph))); + realPhaseOptFourier->SetBinContent(i+1, scale*(re->GetBinContent(i+1)*cos(ph) - im->GetBinContent(i+1)*sin(ph))); realPhaseOptFourier->SetBinError(i+1, 0.0); } diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp index 0b9cbc1f..17e06de2 100644 --- a/src/classes/PFourierCanvas.cpp +++ b/src/classes/PFourierCanvas.cpp @@ -341,6 +341,16 @@ void PFourierCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *select CleanupAverage(); PlotFourier(); } + } else if (x == 'c') { + Int_t state = fFourierPad->GetCrosshair(); + if (state == 0) { + fMainCanvas->ToggleEventStatus(); + fFourierPad->SetCrosshair(2); + } else { + fMainCanvas->ToggleEventStatus(); + fFourierPad->SetCrosshair(0); + } + fMainCanvas->Update(); } else if (x == '+') { // increment phase (Fourier real/imag) IncrementFourierPhase(); } else if (x == '-') { // decrement phase (Fourier real/imag) @@ -366,32 +376,62 @@ void PFourierCanvas::HandleMenuPopup(Int_t id) fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL); fCurrentPlotView = FOURIER_PLOT_REAL; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG); fCurrentPlotView = FOURIER_PLOT_IMAG; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG); fCurrentPlotView = P_MENU_ID_FOURIER_REAL_AND_IMAG; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR); fCurrentPlotView = FOURIER_PLOT_POWER; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE); fCurrentPlotView = FOURIER_PLOT_PHASE; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_OPT_REAL) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_OPT_REAL); fCurrentPlotView = FOURIER_PLOT_PHASE_OPT_REAL; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_PLUS) { IncrementFourierPhase(); } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_MINUS) { @@ -814,7 +854,10 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierIm = fFourier[i]->GetImaginaryFourier(); fFourierHistos[i].dataFourierPwr = fFourier[i]->GetPowerFourier(); fFourierHistos[i].dataFourierPhase = fFourier[i]->GetPhaseFourier(); - fFourierHistos[i].dataFourierPhaseOptReal = fFourier[i]->GetPhaseOptRealFourier(fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + fFourierHistos[i].dataFourierPhaseOptReal = fFourier[i]->GetPhaseOptRealFourier(fFourierHistos[i].dataFourierRe, + fFourierHistos[i].dataFourierIm, fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); + } } // rescale histo to abs(maximum) == 1 @@ -879,17 +922,19 @@ void PFourierCanvas::InitFourierDataSets() } } - // phase opt real - for (UInt_t i=0; iGetBinContent(j); - if (fabs(dval) > max) - max = dval; + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + // phase opt real + for (UInt_t i=0; iGetBinContent(j); + if (fabs(dval) > max) + max = dval; + } } - } - for (UInt_t i=0; iGetNbinsX(); j++) { - fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max)); + for (UInt_t i=0; iGetNbinsX(); j++) { + fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max)); + } } } @@ -903,8 +948,10 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierPwr->SetLineColor(fColorList[i]); fFourierHistos[i].dataFourierPhase->SetMarkerColor(fColorList[i]); fFourierHistos[i].dataFourierPhase->SetLineColor(fColorList[i]); - fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); - fFourierHistos[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + } } // set the marker symbol and size @@ -917,8 +964,10 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierPwr->SetMarkerSize(0.7); fFourierHistos[i].dataFourierPhase->SetMarkerStyle(fMarkerList[i]); fFourierHistos[i].dataFourierPhase->SetMarkerSize(0.7); - fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerStyle(fMarkerList[i]); - fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerSize(0.7); + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerSize(0.7); + } } // initialize average histos @@ -1116,6 +1165,10 @@ void PFourierCanvas::HandleAverage() TString name(""); Double_t dval=0.0; + Bool_t phaseOptRealPresent = false; + if (fFourierHistos[0].dataFourierPhaseOptReal != 0) + phaseOptRealPresent = true; + // check if ALL data shall be averaged if (fAveragedView) { fFourierAverage.resize(1); @@ -1141,10 +1194,12 @@ void PFourierCanvas::HandleAverage() fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmin(), fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmax()); - name = TString(fFourierHistos[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; - fFourierAverage[0].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + if (phaseOptRealPresent) { + name = TString(fFourierHistos[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fFourierAverage[0].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } // real average for (Int_t j=0; jGetNbinsX(); j++) { @@ -1206,20 +1261,22 @@ void PFourierCanvas::HandleAverage() fFourierAverage[0].dataFourierPhase->SetMarkerSize(0.8); fFourierAverage[0].dataFourierPhase->SetMarkerStyle(22); - // phase optimised real average - for (Int_t j=0; jGetNbinsX(); j++) { - dval = 0.0; - for (UInt_t i=0; iGetNbinsX()) - dval += GetInterpolatedValue(fFourierHistos[i].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + if (phaseOptRealPresent) { + // phase optimised real average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (UInt_t i=0; iGetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[i].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + } + fFourierAverage[0].dataFourierPhaseOptReal->SetBinContent(j, dval/fFourierHistos.size()); } - fFourierAverage[0].dataFourierPhaseOptReal->SetBinContent(j, dval/fFourierHistos.size()); + // set marker color, line color, maker size, marker type + fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerColor(kBlack); + fFourierAverage[0].dataFourierPhaseOptReal->SetLineColor(kBlack); + fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerSize(0.8); + fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerStyle(22); } - // set marker color, line color, maker size, marker type - fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerColor(kBlack); - fFourierAverage[0].dataFourierPhaseOptReal->SetLineColor(kBlack); - fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerSize(0.8); - fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerStyle(22); } // check if per data set shall be averaged @@ -1283,10 +1340,12 @@ void PFourierCanvas::HandleAverage() fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmin(), fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmax()); - name = TString(fFourierHistos[start].dataFourierPhaseOptReal->GetTitle()) + "_avg"; - fFourierAverage[i].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + if (phaseOptRealPresent) { + name = TString(fFourierHistos[start].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fFourierAverage[i].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } // real average for (Int_t j=0; jGetNbinsX(); j++) { @@ -1348,24 +1407,61 @@ void PFourierCanvas::HandleAverage() fFourierAverage[i].dataFourierPhase->SetMarkerSize(0.8); fFourierAverage[i].dataFourierPhase->SetMarkerStyle(22); // closed up triangle - // phase optimised real average - for (Int_t j=0; jGetNbinsX(); j++) { - dval = 0.0; - for (Int_t k=start; k<=end; k++) { - if (j < fFourierHistos[k].dataFourierPhaseOptReal->GetNbinsX()) - dval += GetInterpolatedValue(fFourierHistos[k].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + if (phaseOptRealPresent) { + // phase optimised real average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (Int_t k=start; k<=end; k++) { + if (j < fFourierHistos[k].dataFourierPhaseOptReal->GetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[k].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + } + fFourierAverage[i].dataFourierPhaseOptReal->SetBinContent(j, dval/(end-start+1)); } - fFourierAverage[i].dataFourierPhaseOptReal->SetBinContent(j, dval/(end-start+1)); + // set marker color, line color, maker size, marker type + fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); + fFourierAverage[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerSize(0.8); + fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerStyle(22); // closed up triangle } - // set marker color, line color, maker size, marker type - fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); - fFourierAverage[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); - fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerSize(0.8); - fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerStyle(22); // closed up triangle } } } +//-------------------------------------------------------------------------- +// CalcPhaseOptReal (private) +//-------------------------------------------------------------------------- +/** + *

calculate the phase opt. real FT + */ +void PFourierCanvas::CalcPhaseOptReal() +{ + Int_t start = fFourierHistos[0].dataFourierRe->FindFixBin(fInitialXRange[0]); + Int_t end = fFourierHistos[0].dataFourierRe->FindFixBin(fInitialXRange[1]); + + Double_t dval=0.0, max=0.0; + for (UInt_t i=0; iGetPhaseOptRealFourier(fFourierHistos[i].dataFourierRe, + fFourierHistos[i].dataFourierIm, fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); + // normalize it + max = 0.0; + for (Int_t j=start; j<=end; j++) { + dval = fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + for (Int_t j=1; jGetNbinsX(); j++) { + fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max)); + } + // set the marker and line color + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + + // set the marker symbol and size + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerSize(0.7); + } +} + //-------------------------------------------------------------------------- // PlotFourier (private) //-------------------------------------------------------------------------- @@ -1374,6 +1470,14 @@ void PFourierCanvas::HandleAverage() */ void PFourierCanvas::PlotFourier() { + // check if phase opt real Fourier spectra already exists if requested, + // and if not calculate them first + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + if (fFourierHistos[0].dataFourierPhaseOptReal == 0) { // not yet calculated + CalcPhaseOptReal(); + } + } + fFourierPad->cd(); Double_t xmin=0, xmax=0; @@ -1669,6 +1773,11 @@ void PFourierCanvas::PlotAverage() fFourierAverage[0].dataFourierPhase->Draw("p"); break; case FOURIER_PLOT_PHASE_OPT_REAL: + if (fFourierHistos[0].dataFourierPhaseOptReal == 0) { + cout << "debug> need to calculate phase opt. average first ..." << endl; + CalcPhaseOptReal(); + HandleAverage(); + } fFourierAverage[0].dataFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax); ymin = GetMinimum(fFourierAverage[0].dataFourierPhaseOptReal, xmin, xmax); ymax = GetMaximum(fFourierAverage[0].dataFourierPhaseOptReal, xmin, xmax); diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 2cae4d9e..e335fac5 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -182,11 +182,14 @@ PMusrCanvas::PMusrCanvas() * \param wh height (in pixels) of the canvas. * \param batch flag: if set true, the canvas will not be displayed. This is used when just dumping of a * graphical output file is wished. + * \param fourier flag: if set true, the canvas will present the Fourier view. + * \param avg flag: if set true, the canvas will present the averages data/Fourier view. */ 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 avg) : + fStartWithFourier(fourier), fStartWithAvg(avg), + fBatchMode(batch), fPlotNumber(number) { fTimeout = 0; fTimeoutTimer = 0; @@ -235,15 +238,17 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title, * \param colorList pre-defined list of colors * \param batch flag: if set true, the canvas will not be displayed. This is used when just dumping of a * graphical output file is wished. + * \param fourier flag: if set true, the canvas will present the Fourier view. + * \param avg flag: if set true, the canvas will present the averages data/Fourier view. */ 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), - fPlotNumber(number), fFourier(fourierDefault), - fMarkerList(markerList), fColorList(colorList) + const Bool_t batch, const Bool_t fourier, const Bool_t avg) : + fStartWithFourier(fourier), fStartWithAvg(avg), fBatchMode(batch), + fPlotNumber(number), fFourier(fourierDefault), + fMarkerList(markerList), fColorList(colorList) { fTimeout = 0; fTimeoutTimer = 0; @@ -832,6 +837,12 @@ void PMusrCanvas::UpdateDataTheoryPad() HandleFourier(); PlotFourier(); } + + // if fStartWithAvg=true, start with averaged data/Fourier representation + // fStartWithAvg is given at the command line level + if (fStartWithAvg) { + HandleCmdKey(kKeyPress, (Int_t)'a', 0, 0); + } } //-------------------------------------------------------------------------- @@ -1357,6 +1368,13 @@ void PMusrCanvas::HandleMenuPopup(Int_t id) // set appropriate plot view fPreviousPlotView = fCurrentPlotView; fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + // make sure that phase opt. real indeed exists + if (fData[0].dataFourierPhaseOptReal == 0) { + if (fData[0].dataFourierRe == 0) + HandleFourier(); + else + CalcPhaseOptReFT(); + } // uncheck data fPopupMain->UnCheckEntry(P_MENU_ID_DATA+P_MENU_PLOT_OFFSET*fPlotNumber); // check appropriate fourier popup item @@ -3331,10 +3349,7 @@ void PMusrCanvas::HandleDifference() */ void PMusrCanvas::HandleFourier() { - PDoubleVector phaseParam; - Double_t min=-1.0, max=-1.0; - Double_t re, im, ph; - char hName[1024]; + Double_t re, im; // check if plot type is appropriate for fourier if (fPlotType == MSR_PLOT_NON_MUSR) @@ -3347,9 +3362,9 @@ void PMusrCanvas::HandleFourier() double endTime = fXmax; if (!fStartWithFourier) { // fHistoFrame present, hence get start/end from it bin = fHistoFrame->GetXaxis()->GetFirst(); - startTime = fHistoFrame->GetBinCenter(bin); + startTime = fHistoFrame->GetBinLowEdge(bin); bin = fHistoFrame->GetXaxis()->GetLast(); - endTime = fHistoFrame->GetBinCenter(bin); + endTime = fHistoFrame->GetBinLowEdge(bin)+fHistoFrame->GetBinWidth(bin); } for (UInt_t i=0; iGetMsrFourierList()->fRangeForPhaseCorrection[0]; - max = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; - // eventually min/max needs to be extracted from the 'range_for_phase_correction' parameter of the msr-file - fData[i].dataFourierPhaseOptReal = fourierData.GetPhaseOptRealFourier(phaseParam, scale, min, max); - // set marker and line color fData[i].dataFourierRe->SetMarkerColor(fData[i].data->GetMarkerColor()); fData[i].dataFourierRe->SetLineColor(fData[i].data->GetLineColor()); @@ -3385,22 +3394,18 @@ void PMusrCanvas::HandleFourier() fData[i].dataFourierPwr->SetLineColor(fData[i].data->GetLineColor()); fData[i].dataFourierPhase->SetMarkerColor(fData[i].data->GetMarkerColor()); fData[i].dataFourierPhase->SetLineColor(fData[i].data->GetLineColor()); - fData[i].dataFourierPhaseOptReal->SetMarkerColor(fData[i].data->GetMarkerColor()); - fData[i].dataFourierPhaseOptReal->SetLineColor(fData[i].data->GetLineColor()); // set marker size fData[i].dataFourierRe->SetMarkerSize(1); fData[i].dataFourierIm->SetMarkerSize(1); fData[i].dataFourierPwr->SetMarkerSize(1); fData[i].dataFourierPhase->SetMarkerSize(1); - fData[i].dataFourierPhaseOptReal->SetMarkerSize(1); // set marker type fData[i].dataFourierRe->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierIm->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPwr->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPhase->SetMarkerStyle(fData[i].data->GetMarkerStyle()); - fData[i].dataFourierPhaseOptReal->SetMarkerStyle(fData[i].data->GetMarkerStyle()); // calculate fourier transform of the theory Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3; @@ -3418,31 +3423,18 @@ void PMusrCanvas::HandleFourier() // get power part of the data fData[i].theoryFourierPwr = fourierTheory.GetPowerFourier(scale); // get phase part of the data - fData[i].theoryFourierPhase = fourierTheory.GetPhaseFourier(); - - // get phase optimized real fourier from data - - // clone theory Re FT - strcpy(hName, fData[i].theoryFourierPhase->GetName()); - strcat(hName, "_Opt_Real"); - fData[i].theoryFourierPhaseOptReal = (TH1F*) fData[i].theoryFourierRe->Clone(hName); - - // rotate the theory according to the optimized phase parameters - // first find minBin for min of the phase correction - Int_t minBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindBin(min); - Int_t maxBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindBin(max); - for (Int_t j=1; jGetNbinsX(); j++) { - ph = phaseParam[0] + phaseParam[1] * (Double_t)(j-minBin+1) / (Double_t)(maxBin-minBin); - re = fData[i].theoryFourierRe->GetBinContent(j) * cos(ph) - fData[i].theoryFourierIm->GetBinContent(j) * sin(ph); - fData[i].theoryFourierPhaseOptReal->SetBinContent(j, re); - } + fData[i].theoryFourierPhase = fourierTheory.GetPhaseFourier(); // set line colors for the theory fData[i].theoryFourierRe->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierIm->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPwr->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPhase->SetLineColor(fData[i].theory->GetLineColor()); - fData[i].theoryFourierPhaseOptReal->SetLineColor(fData[i].theory->GetLineColor()); + } + + // phase opt. real FT requested initially in the msr-file, hence calculate it here + if (fCurrentPlotView == PV_FOURIER_PHASE_OPT_REAL) { + CalcPhaseOptReFT(); } // apply global phase if present @@ -3475,42 +3467,6 @@ void PMusrCanvas::HandleFourier() } } } - -/* as: will be obsolate ... - // find optimal Fourier phase if range is given - if ((fFourier.fRangeForPhaseCorrection[0] != -1.0) && (fFourier.fRangeForPhaseCorrection[1] != -1.0)) { - - fCurrentFourierPhase = FindOptimalFourierPhase(); - - // apply optimal Fourier phase - double re, im; - const double cp = TMath::Cos(fCurrentFourierPhase/180.0*TMath::Pi()); - const double sp = TMath::Sin(fCurrentFourierPhase/180.0*TMath::Pi()); - - for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set - // calculate new fourier data set value - re = fData[i].dataFourierRe->GetBinContent(j) * cp + fData[i].dataFourierIm->GetBinContent(j) * sp; - im = fData[i].dataFourierIm->GetBinContent(j) * cp - fData[i].dataFourierRe->GetBinContent(j) * sp; - // overwrite fourier data set value - fData[i].dataFourierRe->SetBinContent(j, re); - fData[i].dataFourierIm->SetBinContent(j, im); - } - } - if ((fData[i].theoryFourierRe != 0) && (fData[i].theoryFourierIm != 0)) { - for (Int_t j=0; jGetNbinsX(); j++) { // loop over a fourier data set - // calculate new fourier data set value - re = fData[i].theoryFourierRe->GetBinContent(j) * cp + fData[i].theoryFourierIm->GetBinContent(j) * sp; - im = fData[i].theoryFourierIm->GetBinContent(j) * cp - fData[i].theoryFourierRe->GetBinContent(j) * sp; - // overwrite fourier data set value - fData[i].theoryFourierRe->SetBinContent(j, re); - fData[i].theoryFourierIm->SetBinContent(j, im); - } - } - } - } -*/ } } @@ -4190,6 +4146,10 @@ void PMusrCanvas::CleanupFourierDifference() delete fData[i].diffFourierPhase; fData[i].diffFourierPhase = 0; } + if (fData[i].diffFourierPhaseOptReal != 0) { + delete fData[i].diffFourierPhaseOptReal; + fData[i].diffFourierPhaseOptReal = 0; + } } } @@ -4275,6 +4235,64 @@ void PMusrCanvas::CleanupAverage() } } +//-------------------------------------------------------------------------- +// CalculateDiff (private) +//-------------------------------------------------------------------------- +/** + * @brief PMusrCanvas::CalcPhaseOptReFT + */ +void PMusrCanvas::CalcPhaseOptReFT() +{ + Double_t min = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; + Double_t max = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; + + if ((min == -1.0) && (max == -1.0)) { + if ((fFourier.fPlotRange[0] != -1) && (fFourier.fPlotRange[1] != -1)) { + min = fFourier.fPlotRange[0]; + max = fFourier.fPlotRange[1]; + } else { + min = fData[0].dataFourierRe->GetBinLowEdge(1); + max = fData[0].dataFourierRe->GetBinLowEdge(fData[0].dataFourierRe->GetNbinsX())+fData[0].dataFourierRe->GetBinWidth(1); + } + } + + PDoubleVector phaseParam; + Char_t hName[1024]; + Double_t ph, re; + + for (UInt_t i=0; iSetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPhaseOptReal->SetLineColor(fData[i].data->GetLineColor()); + // set marker size + fData[i].dataFourierPhaseOptReal->SetMarkerSize(1); + // set marker type + fData[i].dataFourierPhaseOptReal->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + + // handle Fourier theory part + // clone theory Re FT + strcpy(hName, fData[i].theoryFourierPhase->GetName()); + strcat(hName, "_Opt_Real"); + fData[i].theoryFourierPhaseOptReal = (TH1F*) fData[i].theoryFourierRe->Clone(hName); + + // rotate the theory according to the optimized phase parameters + // first find minBin for min of the phase correction + Int_t minBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindFixBin(min); + Int_t maxBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindFixBin(max); + + for (Int_t j=1; jGetNbinsX(); j++) { + ph = phaseParam[0] + phaseParam[1] * (Double_t)(j-minBin+1) / (Double_t)(maxBin-minBin); + re = fData[i].theoryFourierRe->GetBinContent(j) * cos(ph) - fData[i].theoryFourierIm->GetBinContent(j) * sin(ph); + fData[i].theoryFourierPhaseOptReal->SetBinContent(j, re); + } + // set line colors for the theory + fData[i].theoryFourierPhaseOptReal->SetLineColor(fData[i].theory->GetLineColor()); + } +} + //-------------------------------------------------------------------------- // CalculateDiff (private) //-------------------------------------------------------------------------- @@ -5924,6 +5942,58 @@ void PMusrCanvas::PlotFourierDifference(Bool_t unzoom) PlotFourierPhaseValue(); + break; + case PV_FOURIER_PHASE_OPT_REAL: + // set x-range + if ((fFourier.fPlotRange[0] != -1) && (fFourier.fPlotRange[1] != -1)) { + xmin = fFourier.fPlotRange[0]; + xmax = fFourier.fPlotRange[1]; + } else { + xmin = fData[0].diffFourierPhaseOptReal->GetBinLowEdge(1); + xmax = fData[0].diffFourierPhaseOptReal->GetBinLowEdge(fData[0].diffFourierPhaseOptReal->GetNbinsX())+fData[0].diffFourierPhaseOptReal->GetBinWidth(1); + } + + // set y-range + // first find minimum/maximum of all histos + ymin = GetMinimum(fData[0].diffFourierPhaseOptReal); + ymax = GetMaximum(fData[0].diffFourierPhaseOptReal); + for (UInt_t i=1; i ymax) + ymax = binContent; + } + + // delete old fHistoFrame if present + if (fHistoFrame) { + delete fHistoFrame; + fHistoFrame = 0; + } + + fHistoFrame = fDataTheoryPad->DrawFrame(xmin, 1.05*ymin, xmax, 1.05*ymax); + + // set ranges for phase opt. real Fourier difference + for (UInt_t i=0; iGetXaxis()->SetRangeUser(xmin, xmax); + fData[i].diffFourierPhaseOptReal->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + } + + // set x-axis title + fHistoFrame->GetXaxis()->SetTitle(xAxisTitle.Data()); + + // set y-axis title + fHistoFrame->GetYaxis()->SetTitleOffset(1.3); + if (fData[0].diffFourierTag == 1) + fHistoFrame->GetYaxis()->SetTitle("Real Fourier (d-f: data-theory)"); + else + fHistoFrame->GetYaxis()->SetTitle("Real Fourier (f-d: [(F data)-(F theory)]"); + + // plot data + for (UInt_t i=0; iDraw("plsame"); + } break; default: break; diff --git a/src/include/PFourier.h b/src/include/PFourier.h index 38f7e7ee..5edaddb9 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -108,11 +108,14 @@ class PFourier virtual Double_t GetResolution() { return fResolution; } virtual Double_t GetMaxFreq(); virtual TH1F* GetRealFourier(const Double_t scale = 1.0); - virtual TH1F* GetPhaseOptRealFourier(vector &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); +//as virtual TH1F* GetPhaseOptRealFourier(vector &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0); virtual TH1F* GetPowerFourier(const Double_t scale = 1.0); virtual TH1F* GetPhaseFourier(const Double_t scale = 1.0); + static TH1F* GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, vector &phase, + const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); + virtual Bool_t IsValid() { return fValid; } private: @@ -136,7 +139,7 @@ class PFourier fftw_complex *fIn; ///< real part of the Fourier transform fftw_complex *fOut; ///< imaginary part of the Fourier transform - PFTPhaseCorrection *fPhCorrectedReFT; +//as PFTPhaseCorrection *fPhCorrectedReFT; virtual void PrepareFFTwInputData(UInt_t apodizationTag); virtual void ApodizeData(Int_t apodizationTag); diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h index 3aff49e9..05dad1e9 100644 --- a/src/include/PFourierCanvas.h +++ b/src/include/PFourierCanvas.h @@ -159,6 +159,7 @@ class PFourierCanvas : public TObject, public TQObject virtual void InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh); virtual void CleanupAverage(); virtual void HandleAverage(); + virtual void CalcPhaseOptReal(); virtual void PlotFourier(); virtual void PlotFourierPhaseValue(); diff --git a/src/musrFT.cpp b/src/musrFT.cpp index 84bc5991..85ca609f 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -713,7 +713,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); } @@ -721,7 +721,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); @@ -1020,7 +1020,6 @@ Int_t main(Int_t argc, Char_t *argv[]) PStartupOptions startup_options; startup_options.writeExpectedChisq = false; startup_options.estimateN0 = true; - startup_options.alphaEstimateN0 = 0.0; TSAXParser *saxParser = new TSAXParser(); PStartupHandler *startupHandler = new PStartupHandler(); if (!startupHandler->StartupFileFound()) { @@ -1057,7 +1056,6 @@ Int_t main(Int_t argc, Char_t *argv[]) } } } - startupHandler->SetStartupOptions(startup_options); // defines the raw time-domain data vector PPrepFourier data(startupParam.packing, startupParam.bkg_range, startupParam.bkg); @@ -1066,7 +1064,7 @@ Int_t main(Int_t argc, Char_t *argv[]) vector msrHandler; msrHandler.resize(startupParam.msrFln.size()); for (UInt_t i=0; iGetStartupOptions(), true); + msrHandler[i] = new PMsrHandler(startupParam.msrFln[i].Data(), &startup_options, true); status = msrHandler[i]->ReadMsrFile(); if (status != PMUSR_SUCCESS) { switch (status) { diff --git a/src/musrview.cpp b/src/musrview.cpp index f91dd6fd..c416b71d 100644 --- a/src/musrview.cpp +++ b/src/musrview.cpp @@ -63,6 +63,7 @@ void musrview_syntax() cout << endl << " --help : display this help and exit."; cout << endl << " --version : output version information and exit."; cout << endl << " -f, --fourier: will directly present the Fourier transform of the ."; + cout << endl << " -a, --avg: will directly present the averaged data/Fourier of the ."; cout << endl << " --: "; cout << endl << " will produce a graphics-output-file without starting a root session."; cout << endl << " the name is based on the , e.g. 3310.msr -> 3310_0.png"; @@ -102,6 +103,7 @@ int main(int argc, char *argv[]) bool success = true; char fileName[128]; bool fourier = false; + bool avg = false; bool graphicsOutput = false; bool asciiOutput = false; char graphicsExtension[128]; @@ -135,6 +137,8 @@ int main(int argc, char *argv[]) break; } else if (!strcmp(argv[i], "-f") || !strcmp(argv[i], "--fourier")) { fourier = true; + } else if (!strcmp(argv[i], "-a") || !strcmp(argv[i], "--avg")) { + avg = true; } else if (!strcmp(argv[i], "--eps") || !strcmp(argv[i], "--pdf") || !strcmp(argv[i], "--gif") || !strcmp(argv[i], "--jpg") || !strcmp(argv[i], "--png") || !strcmp(argv[i], "--svg") || !strcmp(argv[i], "--xpm") || !strcmp(argv[i], "--root")) { @@ -308,12 +312,12 @@ int main(int argc, char *argv[]) startupHandler->GetMarkerList(), startupHandler->GetColorList(), graphicsOutput||asciiOutput, - fourier); + fourier, avg); else musrCanvas = new PMusrCanvas(i, msrHandler->GetMsrTitle()->Data(), 10+i*100, 10+i*100, 800, 600, graphicsOutput||asciiOutput, - fourier); + fourier, avg); if (!musrCanvas->IsValid()) { cerr << endl << ">> musrview **SEVERE ERROR** Couldn't invoke all necessary objects, will quit."; From 5d80c342e0124ad032a98b5635cd19dc0f4bcb10 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 18 Dec 2016 10:39:43 +0100 Subject: [PATCH 10/13] cleaned up some inconsistencies in the musrfit_startup.xml options --- src/classes/PStartupHandler.cpp | 25 ------------------- src/include/PMusr.h | 3 +-- src/include/PStartupHandler.h | 9 +++---- src/musrfit.cpp | 44 ++++----------------------------- src/musrfit_startup.xml | 5 ---- 5 files changed, 9 insertions(+), 77 deletions(-) diff --git a/src/classes/PStartupHandler.cpp b/src/classes/PStartupHandler.cpp index 5660ea88..6ce00364 100644 --- a/src/classes/PStartupHandler.cpp +++ b/src/classes/PStartupHandler.cpp @@ -178,10 +178,6 @@ void PStartupHandler::OnStartDocument() fFourierDefaults.fPlotRange[0] = -1.0; fFourierDefaults.fPlotRange[1] = -1.0; fFourierDefaults.fPhaseIncrement = 1.0; - - fStartupOptions.writeExpectedChisq = false; - fStartupOptions.estimateN0 = true; - fStartupOptions.alphaEstimateN0 = 0.0; } //-------------------------------------------------------------------------- @@ -210,12 +206,6 @@ void PStartupHandler::OnStartElement(const Char_t *str, const TList *attributes) { if (!strcmp(str, "data_path")) { fKey = eDataPath; - } else if (!strcmp(str, "write_per_run_block_chisq")) { - fKey = eWritePerRunBlockChisq; - } else if (!strcmp(str, "estimate_n0")) { - fKey = eEstimateN0; - } else if (!strcmp(str, "alpha_estimate_n0")) { - fKey = eAlphaEstimateN0; } else if (!strcmp(str, "marker")) { fKey = eMarker; } else if (!strcmp(str, "color")) { @@ -270,21 +260,6 @@ void PStartupHandler::OnCharacters(const Char_t *str) // add str to the path list fDataPathList.push_back(str); break; - case eWritePerRunBlockChisq: - tstr = TString(str); - if (tstr.BeginsWith("y") || tstr.BeginsWith("Y")) - fStartupOptions.writeExpectedChisq = true; - break; - case eEstimateN0: - tstr = TString(str); - if (tstr.BeginsWith("n") || tstr.BeginsWith("N")) - fStartupOptions.estimateN0 = false; - break; - case eAlphaEstimateN0: - tstr = TString(str); - if (tstr.IsFloat()) - fStartupOptions.alphaEstimateN0 = tstr.Atof(); - 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 96df487f..d093cbfa 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -815,12 +815,11 @@ typedef struct { //------------------------------------------------------------- /** - *

Holds the informations for the any2many converter program + *

Holds the informations */ 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 } PStartupOptions; //------------------------------------------------------------- diff --git a/src/include/PStartupHandler.h b/src/include/PStartupHandler.h index 9125a67b..89602014 100644 --- a/src/include/PStartupHandler.h +++ b/src/include/PStartupHandler.h @@ -73,17 +73,15 @@ class PStartupHandler : public TObject, public TQObject virtual void CheckLists(); - virtual PStartupOptions* GetStartupOptions() { return &fStartupOptions; } ///< returns the startup options virtual PMsrFourierStructure GetFourierDefaults() { return fFourierDefaults; } ///< returns the Fourier defaults virtual const PStringVector GetDataPathList() const { return fDataPathList; } ///< returns the search data path list virtual const PIntVector GetMarkerList() const { return fMarkerList; } ///< returns the marker list virtual const PIntVector GetColorList() const { return fColorList; } ///< returns the color list - virtual void SetStartupOptions(const PStartupOptions opt) { fStartupOptions = opt; } - private: - enum EKeyWords {eEmpty, eComment, eDataPath, eOptions, eWritePerRunBlockChisq, eEstimateN0, eAlphaEstimateN0, - eFourierSettings, eUnits, eFourierPower, eApodization, ePlot, ePhase, ePhaseIncrement, + enum EKeyWords {eEmpty, eComment, eDataPath, eOptions, + eFourierSettings, eUnits, eFourierPower, + eApodization, ePlot, ePhase, ePhaseIncrement, eRootSettings, eMarkerList, eMarker, eColorList, eColor}; EKeyWords fKey; ///< xml filter key @@ -94,7 +92,6 @@ class PStartupHandler : public TObject, public TQObject PStringVector fDataPathList; ///< search data path list PIntVector fMarkerList; ///< marker list PIntVector fColorList; ///< color list - PStartupOptions fStartupOptions; ///< collects all startup options which will be requested by PMsrFileHandler Bool_t StartupFileExists(Char_t *fln); diff --git a/src/musrfit.cpp b/src/musrfit.cpp index 0bb0bd2d..aa6200a1 100644 --- a/src/musrfit.cpp +++ b/src/musrfit.cpp @@ -104,9 +104,8 @@ void musrfit_syntax() cout << endl << " -t, --title-from-data-file: will replace the run title by the"; cout << endl << " run title of the FIRST run of the run block, if a run title"; cout << endl << " is present in the data file."; - cout << endl << " -e, --estimateN0: estimate N0 for single histogram fits flag. can have the values 'yes' or 'no'."; + cout << endl << " -e, --estimateN0: estimate N0 for single histogram fits."; cout << endl << " -p, --per-run-block-chisq: will write per run block chisq to the msr-file."; - cout << endl << " can have the values 'yes' or 'no'."; cout << endl << " --dump is writing a data file with the fit data and the theory"; cout << endl << " can be 'ascii', 'root'"; cout << endl << " --timeout : overwrites to predefined timeout of " << timeout << " (sec)."; @@ -444,8 +443,7 @@ int main(int argc, char *argv[]) bool timeout_enabled = true; PStartupOptions startup_options; startup_options.writeExpectedChisq = false; - startup_options.estimateN0 = true; - startup_options.alphaEstimateN0 = 0.0; + startup_options.estimateN0 = false; TString dump(""); char filename[1024]; @@ -496,39 +494,9 @@ int main(int argc, char *argv[]) break; } } else if (!strcmp(argv[i], "-e") || !strcmp(argv[i], "--estimateN0")) { - if (i with unsupported = " << argv[i+1] << endl; - show_syntax = true; - break; - } - i++; - } else { - cerr << endl << "musrfit: **ERROR** found option --estimateN0 without " << endl; - show_syntax = true; - break; - } + startup_options.estimateN0 = true; } else if (!strcmp(argv[i], "-p") || !strcmp(argv[i], "--per-run-block-chisq")) { - if (i with unsupported = " << argv[i+1] << endl; - show_syntax = true; - break; - } - i++; - } else { - cerr << endl << "musrfit: **ERROR** found option --per-run-block-chisq without " << endl; - show_syntax = true; - break; - } + startup_options.writeExpectedChisq = true; } else if (!strcmp(argv[i], "--timeout")) { if (iSetStartupOptions(startup_options); // read msr-file PMsrHandler *msrHandler = 0; if (startupHandler) - msrHandler = new PMsrHandler(filename, startupHandler->GetStartupOptions()); + msrHandler = new PMsrHandler(filename, &startup_options); else msrHandler = new PMsrHandler(filename); status = msrHandler->ReadMsrFile(); diff --git a/src/musrfit_startup.xml b/src/musrfit_startup.xml index daf13703..3aa5ab57 100644 --- a/src/musrfit_startup.xml +++ b/src/musrfit_startup.xml @@ -14,11 +14,6 @@ /afs/psi.ch/project/bulkmusr/data/alc /afs/psi.ch/project/bulkmusr/data/hifi /afs/psi.ch/project/bulkmusr/data/lem - - n - yes - 0.0 - Gauss 0 From 265fe18f36d1d809649f23c562eb354d5d9195e1 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 18 Dec 2016 10:41:34 +0100 Subject: [PATCH 11/13] adopted to the new features in musrview and some additional little improvements. --- src/musredit/PAdmin.cpp | 16 +++++++++ src/musredit/PAdmin.h | 6 +++- src/musredit/PPrefsDialog.cpp | 1 + src/musredit/PPrefsDialog.h | 1 + src/musredit/PTextEdit.cpp | 20 ++++++------ src/musredit/forms/PPrefsDialog.ui | 15 ++++++++- src/musredit_qt5/PAdmin.cpp | 45 ++++++++++++++++++++++---- src/musredit_qt5/PAdmin.h | 6 +++- src/musredit_qt5/PPrefsDialog.cpp | 1 + src/musredit_qt5/PPrefsDialog.h | 1 + src/musredit_qt5/PTextEdit.cpp | 36 +++++++++++++++------ src/musredit_qt5/forms/PPrefsDialog.ui | 15 ++++++++- 12 files changed, 133 insertions(+), 30 deletions(-) diff --git a/src/musredit/PAdmin.cpp b/src/musredit/PAdmin.cpp index d4368c2c..2ab1e9b2 100644 --- a/src/musredit/PAdmin.cpp +++ b/src/musredit/PAdmin.cpp @@ -87,6 +87,8 @@ bool PAdminXMLParser::startElement( const QString&, const QString&, fKeyWord = eTitleFromDataFile; } else if (qName == "musrview_show_fourier") { fKeyWord = eMusrviewShowFourier; + } else if (qName == "musrview_show_avg") { + fKeyWord = eMusrviewShowAvg; } else if (qName == "enable_musrt0") { fKeyWord = eEnableMusrT0; } else if (qName == "keep_minuit2_output") { @@ -251,6 +253,13 @@ bool PAdminXMLParser::characters(const QString& str) flag = false; fAdmin->setMusrviewShowFourierFlag(flag); break; + case eMusrviewShowAvg: + if (str == "y") + flag = true; + else + flag = false; + fAdmin->setMusrviewShowAvgFlag(flag); + break; case eEnableMusrT0: if (str == "y") flag = true; @@ -614,6 +623,7 @@ PAdmin::PAdmin() : QObject() fFileFormat = QString(""); fMusrviewShowFourier = false; + fMusrviewShowAvg = false; fTitleFromDataFile = false; fEnableMusrT0 = false; @@ -837,6 +847,12 @@ int PAdmin::savePrefs(QString pref_fln) else data[i] = " n"; } + if (data[i].contains("") && data[i].contains("")) { + if (fMusrviewShowAvg) + data[i] = " y"; + else + data[i] = " n"; + } if (data[i].contains("") && data[i].contains("")) { if (fEnableMusrT0) data[i] = " y"; diff --git a/src/musredit/PAdmin.h b/src/musredit/PAdmin.h index 861dde15..2d0dc51d 100644 --- a/src/musredit/PAdmin.h +++ b/src/musredit/PAdmin.h @@ -69,7 +69,8 @@ class PAdminXMLParser : public QXmlDefaultHandler private: enum EAdminKeyWords {eEmpty, eTimeout, eKeepMinuit2Output, eDumpAscii, eDumpRoot, - eTitleFromDataFile, eChisqPreRunBlock, eEstimateN0, eMusrviewShowFourier, eEnableMusrT0, + eTitleFromDataFile, eChisqPreRunBlock, eEstimateN0, + eMusrviewShowFourier, eMusrviewShowAvg, eEnableMusrT0, eFontName, eFontSize, eExecPath, eDefaultSavePath, eRecentFile, eBeamline, eInstitute, eFileFormat, eLifetimeCorrection, eMsrDefaultFilePath, eTheoFuncPixmapPath, eFunc, eFuncName, eFuncComment, eFuncLabel, @@ -119,6 +120,7 @@ class PAdmin : public QObject QString getDefaultSavePath() { return fDefaultSavePath; } bool getTitleFromDataFileFlag() { return fTitleFromDataFile; } bool getMusrviewShowFourierFlag() { return fMusrviewShowFourier; } + bool getMusrviewShowAvgFlag() { return fMusrviewShowAvg; } bool getEnableMusrT0Flag() { return fEnableMusrT0; } bool getKeepMinuit2OutputFlag() { return fKeepMinuit2Output; } bool getDumpAsciiFlag() { return fDumpAscii; } @@ -142,6 +144,7 @@ class PAdmin : public QObject void setTimeout(const int ival) { fTimeout = ival; } void setTitleFromDataFileFlag(const bool flag) { fTitleFromDataFile = flag; } void setMusrviewShowFourierFlag(const bool flag) { fMusrviewShowFourier = flag; } + void setMusrviewShowAvgFlag(const bool flag) { fMusrviewShowAvg = flag; } void setEnableMusrT0Flag(const bool flag) { fEnableMusrT0 = flag; } void setKeepMinuit2OutputFlag(const bool flag) { fKeepMinuit2Output = flag; } void setDumpAsciiFlag(const bool flag) { fDumpAscii = flag; } @@ -185,6 +188,7 @@ class PAdmin : public QObject QVector fRecentFile; ///< keep vector of recent path-file names bool fMusrviewShowFourier; ///< flag indicating if musrview should show at startup data (=false) or Fourier of data (=true). + bool fMusrviewShowAvg; ///< flag indicating if musrview should show at startup averaged (=true) or original (=false) data/Fourier. bool fKeepMinuit2Output; ///< flag indicating if the Minuit2 output shall be kept (default: no) bool fDumpAscii; ///< flag indicating if musrfit shall make an ascii-dump file (for debugging purposes, default: no). bool fDumpRoot; ///< flag indicating if musrfit shall make an root-dump file (for debugging purposes, default: no). diff --git a/src/musredit/PPrefsDialog.cpp b/src/musredit/PPrefsDialog.cpp index c876a12f..1b40a423 100644 --- a/src/musredit/PPrefsDialog.cpp +++ b/src/musredit/PPrefsDialog.cpp @@ -62,6 +62,7 @@ PPrefsDialog::PPrefsDialog(PAdmin *admin) : fAdmin(admin) fPerRunBlockChisq_checkBox->setChecked(fAdmin->getChisqPerRunBlockFlag()); fEstimateN0_checkBox->setChecked(fAdmin->getEstimateN0Flag()); fFourier_checkBox->setChecked(fAdmin->getMusrviewShowFourierFlag()); + fAvg_checkBox->setChecked(fAdmin->getMusrviewShowAvgFlag()); fTimeout_lineEdit->setText(QString("%1").arg(fAdmin->getTimeout())); fTimeout_lineEdit->setValidator(new QIntValidator(fTimeout_lineEdit)); diff --git a/src/musredit/PPrefsDialog.h b/src/musredit/PPrefsDialog.h index 7543c0a8..edba329f 100644 --- a/src/musredit/PPrefsDialog.h +++ b/src/musredit/PPrefsDialog.h @@ -47,6 +47,7 @@ class PPrefsDialog : public QDialog, private Ui::PPrefsDialog PPrefsDialog(PAdmin *admin); bool getMusrviewShowFourierFlag() { return fFourier_checkBox->isChecked(); } + bool getMusrviewShowAvgFlag() { return fAvg_checkBox->isChecked(); } bool getKeepMinuit2OutputFlag() { return fKeepMn2Output_checkBox->isChecked(); } bool getTitleFromDataFileFlag() { return fTitleFromData_checkBox->isChecked(); } bool getEnableMusrT0Flag() { return fEnableMusrT0_checkBox->isChecked(); } diff --git a/src/musredit/PTextEdit.cpp b/src/musredit/PTextEdit.cpp index 5774fa8b..6bf78ca2 100644 --- a/src/musredit/PTextEdit.cpp +++ b/src/musredit/PTextEdit.cpp @@ -1699,6 +1699,12 @@ void PTextEdit::musrCalcChisq() if ( !currentEditor() ) return; + int result = 0; + if (fAdmin->getEstimateN0Flag()) + result = QMessageBox::question(this, "Estimate N0 active", + "Do you wish a chisq/mlh evaluation with an automatic N0 estimate?", + QMessageBox::Yes, QMessageBox::No); + QString tabLabel = fTabWidget->tabText(fTabWidget->currentIndex()); if (tabLabel == "noname") { QMessageBox::critical(this, "**ERROR**", "For a fit a real msr-file is needed."); @@ -1716,8 +1722,8 @@ void PTextEdit::musrCalcChisq() cmd.append(str); cmd.append(QFileInfo(*fFilenames.find( currentEditor())).fileName() ); cmd.append("--chisq-only"); - cmd.append("--estimateN0"); - cmd.append("no"); + if (fAdmin->getEstimateN0Flag() && (result == QMessageBox::Yes)) + cmd.append("--estimateN0"); PFitOutputHandler fitOutputHandler(QFileInfo(*fFilenames.find( currentEditor() )).absolutePath(), cmd); fitOutputHandler.setModal(true); fitOutputHandler.exec(); @@ -1770,19 +1776,11 @@ void PTextEdit::musrFit() // check estimate N0 flag if (fAdmin->getEstimateN0Flag()) { cmd.append("--estimateN0"); - cmd.append("yes"); - } else { - cmd.append("--estimateN0"); - cmd.append("no"); } // check per-run-block-chisq flag if (fAdmin->getChisqPerRunBlockFlag()) { cmd.append("--per-run-block-chisq"); - cmd.append("yes"); - } else { - cmd.append("--per-run-block-chisq"); - cmd.append("no"); } // add timeout @@ -2203,6 +2201,8 @@ void PTextEdit::musrView() cmd += str + "\" --timeout " + numStr; if (fAdmin->getMusrviewShowFourierFlag()) cmd += " -f "; + if (fAdmin->getMusrviewShowAvgFlag()) + cmd += " -a "; cmd += " &"; int status=system(cmd.toLatin1()); diff --git a/src/musredit/forms/PPrefsDialog.ui b/src/musredit/forms/PPrefsDialog.ui index 90f3e0b9..0bfc0f51 100644 --- a/src/musredit/forms/PPrefsDialog.ui +++ b/src/musredit/forms/PPrefsDialog.ui @@ -33,7 +33,7 @@ - 0 + 2 @@ -175,6 +175,19 @@ start with Fourier + + + + 10 + 40 + 261 + 22 + + + + start with averaged data/Fourier + + diff --git a/src/musredit_qt5/PAdmin.cpp b/src/musredit_qt5/PAdmin.cpp index 2f6b1bc4..4a1f2fa1 100644 --- a/src/musredit_qt5/PAdmin.cpp +++ b/src/musredit_qt5/PAdmin.cpp @@ -37,6 +37,8 @@ using namespace std; #include #include +#include + #include "PAdmin.h" //-------------------------------------------------------------------------- @@ -87,6 +89,8 @@ bool PAdminXMLParser::startElement( const QString&, const QString&, fKeyWord = eTitleFromDataFile; } else if (qName == "musrview_show_fourier") { fKeyWord = eMusrviewShowFourier; + } else if (qName == "musrview_show_avg") { + fKeyWord = eMusrviewShowAvg; } else if (qName == "enable_musrt0") { fKeyWord = eEnableMusrT0; } else if (qName == "keep_minuit2_output") { @@ -234,7 +238,7 @@ bool PAdminXMLParser::characters(const QString& str) case eExecPath: fAdmin->setExecPath(QString(str.toLatin1()).trimmed()); break; - case eDefaultSavePath: + case eDefaultSavePath: fAdmin->setDefaultSavePath(QString(str.toLatin1()).trimmed()); break; case eTitleFromDataFile: @@ -251,6 +255,13 @@ bool PAdminXMLParser::characters(const QString& str) flag = false; fAdmin->setMusrviewShowFourierFlag(flag); break; + case eMusrviewShowAvg: + if (str == "y") + flag = true; + else + flag = false; + fAdmin->setMusrviewShowAvgFlag(flag); + break; case eEnableMusrT0: if (str == "y") flag = true; @@ -477,7 +488,7 @@ bool PAdminXMLParser::endDocument() str = expandPath(fAdmin->getDefaultSavePath()); if (!str.isEmpty()) fAdmin->setDefaultSavePath(str); - } + } if (fAdmin->getMsrDefaultFilePath().indexOf('$') >= 0) { str = expandPath(fAdmin->getMsrDefaultFilePath()); @@ -567,13 +578,21 @@ QString PAdminXMLParser::expandPath(const QString &str) QString msg; QString newStr=""; + QProcessEnvironment procEnv = QProcessEnvironment::systemEnvironment(); + QStringList list = str.split("/"); for ( QStringList::Iterator it = list.begin(); it != list.end(); ++it ) { token = *it; if (token.contains("$")) { token.remove('$'); - path = std::getenv(token.toLatin1()); + if (!procEnv.contains(token)) { + msg = QString("Couldn't find '%1'. Some things might not work properly").arg(token); + QMessageBox::warning(0, "**WARNING**", msg, QMessageBox::Ok, QMessageBox::NoButton); + newStr = ""; + break; + } + path = procEnv.value(token, ""); if (path.isEmpty()) { msg = QString("Couldn't expand '%1'. Some things might not work properly").arg(token); QMessageBox::warning(0, "**WARNING**", msg, QMessageBox::Ok, QMessageBox::NoButton); @@ -614,6 +633,7 @@ PAdmin::PAdmin() : QObject() fFileFormat = QString(""); fMusrviewShowFourier = false; + fMusrviewShowAvg = false; fTitleFromDataFile = false; fEnableMusrT0 = false; @@ -646,17 +666,18 @@ PAdmin::PAdmin() : QObject() QString path = QString("./"); QString fln = QString("musredit_startup.xml"); QString pathFln = path + fln; + QProcessEnvironment procEnv = QProcessEnvironment::systemEnvironment(); if (!QFile::exists(pathFln)) { // 2nd: check $HOME/.musrfit/musredit/musredit_startup.xml - path = std::getenv("HOME"); + path = procEnv.value("HOME", ""); pathFln = path + "/.musrfit/musredit/" + fln; if (!QFile::exists(pathFln)) { // 3rd: check $MUSRFITPATH/musredit_startup.xml - path = std::getenv("MUSRFITPATH"); + path = procEnv.value("MUSRFITPATH", ""); pathFln = path + "/" + fln; if (!QFile::exists(pathFln)) { // 4th: check $ROOTSYS/bin/musredit_startup.xml - path = std::getenv("ROOTSYS"); + path = procEnv.value("ROOTSYS", ""); pathFln = path + "/bin/" + fln; } } @@ -837,12 +858,24 @@ int PAdmin::savePrefs(QString pref_fln) else data[i] = " n"; } + if (data[i].contains("") && data[i].contains("")) { + if (fMusrviewShowAvg) + data[i] = " y"; + else + data[i] = " n"; + } if (data[i].contains("") && data[i].contains("")) { if (fEnableMusrT0) data[i] = " y"; else data[i] = " n"; } + if (data[i].contains("") && data[i].contains("")) { + data[i] = QString(" %1").arg(fFontName); + } + if (data[i].contains("") && data[i].contains("")) { + data[i] = QString(" %1").arg(fFontSize); + } } // write prefs diff --git a/src/musredit_qt5/PAdmin.h b/src/musredit_qt5/PAdmin.h index ef9abfc3..b966ef78 100644 --- a/src/musredit_qt5/PAdmin.h +++ b/src/musredit_qt5/PAdmin.h @@ -69,7 +69,8 @@ class PAdminXMLParser : public QXmlDefaultHandler private: enum EAdminKeyWords {eEmpty, eTimeout, eKeepMinuit2Output, eDumpAscii, eDumpRoot, - eTitleFromDataFile, eChisqPreRunBlock, eEstimateN0, eMusrviewShowFourier, eEnableMusrT0, + eTitleFromDataFile, eChisqPreRunBlock, eEstimateN0, + eMusrviewShowFourier, eMusrviewShowAvg, eEnableMusrT0, eFontName, eFontSize, eExecPath, eDefaultSavePath, eRecentFile, eBeamline, eInstitute, eFileFormat, eLifetimeCorrection, eMsrDefaultFilePath, eTheoFuncPixmapPath, eFunc, eFuncName, eFuncComment, eFuncLabel, @@ -119,6 +120,7 @@ class PAdmin : public QObject QString getDefaultSavePath() { return fDefaultSavePath; } bool getTitleFromDataFileFlag() { return fTitleFromDataFile; } bool getMusrviewShowFourierFlag() { return fMusrviewShowFourier; } + bool getMusrviewShowAvgFlag() { return fMusrviewShowAvg; } bool getEnableMusrT0Flag() { return fEnableMusrT0; } bool getKeepMinuit2OutputFlag() { return fKeepMinuit2Output; } bool getDumpAsciiFlag() { return fDumpAscii; } @@ -142,6 +144,7 @@ class PAdmin : public QObject void setTimeout(const int ival) { fTimeout = ival; } void setTitleFromDataFileFlag(const bool flag) { fTitleFromDataFile = flag; } void setMusrviewShowFourierFlag(const bool flag) { fMusrviewShowFourier = flag; } + void setMusrviewShowAvgFlag(const bool flag) { fMusrviewShowAvg = flag; } void setEnableMusrT0Flag(const bool flag) { fEnableMusrT0 = flag; } void setKeepMinuit2OutputFlag(const bool flag) { fKeepMinuit2Output = flag; } void setDumpAsciiFlag(const bool flag) { fDumpAscii = flag; } @@ -185,6 +188,7 @@ class PAdmin : public QObject QVector fRecentFile; ///< keep vector of recent path-file names bool fMusrviewShowFourier; ///< flag indicating if musrview should show at startup data (=false) or Fourier of data (=true). + bool fMusrviewShowAvg; ///< flag indicating if musrview should show at startup averaged (=true) or original (=false) data/Fourier. bool fKeepMinuit2Output; ///< flag indicating if the Minuit2 output shall be kept (default: no) bool fDumpAscii; ///< flag indicating if musrfit shall make an ascii-dump file (for debugging purposes, default: no). bool fDumpRoot; ///< flag indicating if musrfit shall make an root-dump file (for debugging purposes, default: no). diff --git a/src/musredit_qt5/PPrefsDialog.cpp b/src/musredit_qt5/PPrefsDialog.cpp index c876a12f..1b40a423 100644 --- a/src/musredit_qt5/PPrefsDialog.cpp +++ b/src/musredit_qt5/PPrefsDialog.cpp @@ -62,6 +62,7 @@ PPrefsDialog::PPrefsDialog(PAdmin *admin) : fAdmin(admin) fPerRunBlockChisq_checkBox->setChecked(fAdmin->getChisqPerRunBlockFlag()); fEstimateN0_checkBox->setChecked(fAdmin->getEstimateN0Flag()); fFourier_checkBox->setChecked(fAdmin->getMusrviewShowFourierFlag()); + fAvg_checkBox->setChecked(fAdmin->getMusrviewShowAvgFlag()); fTimeout_lineEdit->setText(QString("%1").arg(fAdmin->getTimeout())); fTimeout_lineEdit->setValidator(new QIntValidator(fTimeout_lineEdit)); diff --git a/src/musredit_qt5/PPrefsDialog.h b/src/musredit_qt5/PPrefsDialog.h index 7543c0a8..edba329f 100644 --- a/src/musredit_qt5/PPrefsDialog.h +++ b/src/musredit_qt5/PPrefsDialog.h @@ -47,6 +47,7 @@ class PPrefsDialog : public QDialog, private Ui::PPrefsDialog PPrefsDialog(PAdmin *admin); bool getMusrviewShowFourierFlag() { return fFourier_checkBox->isChecked(); } + bool getMusrviewShowAvgFlag() { return fAvg_checkBox->isChecked(); } bool getKeepMinuit2OutputFlag() { return fKeepMn2Output_checkBox->isChecked(); } bool getTitleFromDataFileFlag() { return fTitleFromData_checkBox->isChecked(); } bool getEnableMusrT0Flag() { return fEnableMusrT0_checkBox->isChecked(); } diff --git a/src/musredit_qt5/PTextEdit.cpp b/src/musredit_qt5/PTextEdit.cpp index abf34026..c0285d1d 100644 --- a/src/musredit_qt5/PTextEdit.cpp +++ b/src/musredit_qt5/PTextEdit.cpp @@ -848,6 +848,13 @@ void PTextEdit::fileOpen() ++it; } + + // in case there is a 1st empty tab "noname", remove it + if (fTabWidget->tabText(0) == "noname") { // has to be the first, otherwise do nothing + fFileSystemWatcher->removePath("noname"); + + delete fTabWidget->widget(0); + } } //---------------------------------------------------------------------------------------------------- @@ -879,6 +886,13 @@ void PTextEdit::fileOpenRecent() else fileReload(); } + + // in case there is a 1st empty tab "noname", remove it + if (fTabWidget->tabText(0) == "noname") { // has to be the first, otherwise do nothing + fFileSystemWatcher->removePath("noname"); + + delete fTabWidget->widget(0); + } } //---------------------------------------------------------------------------------------------------- @@ -1699,6 +1713,11 @@ void PTextEdit::musrCalcChisq() if ( !currentEditor() ) return; + int result = 0; + if (fAdmin->getEstimateN0Flag()) + result = QMessageBox::question(this, "Estimate N0 active", + "Do you wish a chisq/mlh evaluation with an automatic N0 estimate?"); + QString tabLabel = fTabWidget->tabText(fTabWidget->currentIndex()); if (tabLabel == "noname") { QMessageBox::critical(this, "**ERROR**", "For a fit a real msr-file is needed."); @@ -1716,8 +1735,8 @@ void PTextEdit::musrCalcChisq() cmd.append(str); cmd.append(QFileInfo(*fFilenames.find( currentEditor())).fileName() ); cmd.append("--chisq-only"); - cmd.append("--estimateN0"); - cmd.append("no"); + if (fAdmin->getEstimateN0Flag() && (result == QMessageBox::Yes)) + cmd.append("--estimateN0"); PFitOutputHandler fitOutputHandler(QFileInfo(*fFilenames.find( currentEditor() )).absolutePath(), cmd); fitOutputHandler.setModal(true); fitOutputHandler.exec(); @@ -1770,19 +1789,11 @@ void PTextEdit::musrFit() // check estimate N0 flag if (fAdmin->getEstimateN0Flag()) { cmd.append("--estimateN0"); - cmd.append("yes"); - } else { - cmd.append("--estimateN0"); - cmd.append("no"); } // check per-run-block-chisq flag if (fAdmin->getChisqPerRunBlockFlag()) { cmd.append("--per-run-block-chisq"); - cmd.append("yes"); - } else { - cmd.append("--per-run-block-chisq"); - cmd.append("no"); } // add timeout @@ -2205,6 +2216,10 @@ void PTextEdit::musrView() if (fAdmin->getMusrviewShowFourierFlag()) arg << "-f"; + // start with averaged data/Fourier? + if (fAdmin->getMusrviewShowAvgFlag()) + arg << "-a"; + QProcess *proc = new QProcess(this); // make sure that the system environment variables are properly set @@ -2333,6 +2348,7 @@ void PTextEdit::musrPrefs() if (dlg->exec() == QDialog::Accepted) { fAdmin->setMusrviewShowFourierFlag(dlg->getMusrviewShowFourierFlag()); + fAdmin->setMusrviewShowAvgFlag(dlg->getMusrviewShowAvgFlag()); fAdmin->setKeepMinuit2OutputFlag(dlg->getKeepMinuit2OutputFlag()); fAdmin->setTitleFromDataFileFlag(dlg->getTitleFromDataFileFlag()); fAdmin->setEnableMusrT0Flag(dlg->getEnableMusrT0Flag()); diff --git a/src/musredit_qt5/forms/PPrefsDialog.ui b/src/musredit_qt5/forms/PPrefsDialog.ui index 90f3e0b9..8f78eaa4 100644 --- a/src/musredit_qt5/forms/PPrefsDialog.ui +++ b/src/musredit_qt5/forms/PPrefsDialog.ui @@ -17,7 +17,7 @@ Preferences - + :/images/musrfit.xpm:/images/musrfit.xpm @@ -175,6 +175,19 @@ start with Fourier + + + + 10 + 40 + 231 + 21 + + + + start with averaged data/Fourier + + From 43a89451092c18eb2183509a29a52a92cc413ead Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 18 Dec 2016 10:43:00 +0100 Subject: [PATCH 12/13] update of the docu. --- doc/html/user/MUSR/BmwLibs.html | 12 +-- doc/html/user/MUSR/CiteMusrFit.html | 12 +-- doc/html/user/MUSR/LibFitPofB.html | 18 ++--- doc/html/user/MUSR/LibZFRelaxation.html | 22 +++--- doc/html/user/MUSR/Msr2Data.html | 26 +++---- doc/html/user/MUSR/MusrFit.html | 78 +++++++++++++------ .../user/MUSR/MusrFitAcknowledgements.html | 16 ++-- doc/html/user/MUSR/MusrFitSetup.html | 34 ++++---- doc/html/user/MUSR/MusrGui.html | 20 ++--- doc/html/user/MUSR/MusrRoot.html | 26 +++---- doc/html/user/MUSR/QuickStart.html | 20 ++--- doc/html/user/MUSR/TutorialSingleHisto.html | 24 +++--- doc/html/user/MUSR/WebHome.html | 22 +++--- 13 files changed, 182 insertions(+), 148 deletions(-) diff --git a/doc/html/user/MUSR/BmwLibs.html b/doc/html/user/MUSR/BmwLibs.html index 209d2779..051059ca 100644 --- a/doc/html/user/MUSR/BmwLibs.html +++ b/doc/html/user/MUSR/BmwLibs.html @@ -1,6 +1,6 @@ - + @@ -14,14 +14,14 @@ - + - - + - + + - - - - + + + - + - + + + - + + - - - + + + + + - - + + - +