diff --git a/ChangeLog b/ChangeLog
index ac0e3904..7576b4bd 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -4,8 +4,15 @@
changes since 0.14.0
===================================
+NEW 2015-09-24 adding a phase optimized real Fourier to musrFT. This is still VERY
+ experimental, and hence only available from within musrFT. Eventually
+ it should make its way into musrview as well. Furthermore the
+ documentation needs to be written yet. There is only a short not in
+ the command line present.
NEW 2015-02-23 implemented an average-per-data-set option for musrFT.
NEW 2015-02-21 add proper Mac icon to musredit
+FIXED 2015-09-22 generating global fit msr-files including a GLOBAL block is working
+ now as expected.
FIXED 2015-09-17 in PMsr2Data::PrepareGlobalInputFile() there seem to be 'unmotivated'
break; commands in some loops. They prevent a proper map handling.
Since this is a real puzzle I contacted BMW for clarification.
diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp
index 38f94d12..ea08370b 100644
--- a/src/classes/PFourier.cpp
+++ b/src/classes/PFourier.cpp
@@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
- * Copyright (C) 2007-2014 by Andreas Suter *
+ * Copyright (C) 2007-2015 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -163,7 +163,7 @@ PFourier::~PFourier()
}
//--------------------------------------------------------------------------
-// Transform
+// Transform (public)
//--------------------------------------------------------------------------
/**
* Carries out the Fourier transform. It is assumed that fStartTime is the time zero
@@ -201,7 +201,7 @@ void PFourier::Transform(UInt_t apodizationTag)
}
//--------------------------------------------------------------------------
-// GetMaxFreq
+// GetMaxFreq (public)
//--------------------------------------------------------------------------
/**
*
returns the maximal frequency in units choosen, i.e. Gauss, Tesla, MHz, Mc/s
@@ -218,7 +218,7 @@ Double_t PFourier::GetMaxFreq()
}
//--------------------------------------------------------------------------
-// GetRealFourier
+// GetRealFourier (public)
//--------------------------------------------------------------------------
/**
*
returns the real part Fourier as a histogram.
@@ -259,7 +259,108 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
-// GetImaginaryFourier
+// 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$.
+ *
+ * \return the TH1F histo of the phase 'optimzed' real Fourier transform.
+ *
+ * \param phase return value of the optimal phase
+ * \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)
+{
+ UInt_t noOfFourierBins = 0;
+ if (fNoOfBins % 2 == 0)
+ noOfFourierBins = fNoOfBins/2;
+ else
+ noOfFourierBins = (fNoOfBins+1)/2;
+
+ UInt_t minBin = 0;
+ UInt_t maxBin = noOfFourierBins;
+
+ // check if minimum frequency is given. If yes, get the proper minBin
+ if (min != -1.0) {
+ minBin = (UInt_t)(min/fResolution);
+ }
+
+ // check if maximum frequency is given. If yes, get the proper maxBin
+ if (max != -1.0) {
+ maxBin = (UInt_t)(max/fResolution);
+ }
+
+ // 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]);
+ }
+
+ // 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; jGetName());
+ snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", fData->GetTitle());
+
+ TH1F *realPhaseOptFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/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
+ for (UInt_t i=0; iSetBinContent(i+1, scale*(fOut[i][0]*cos(phase) + fOut[i][1]*sin(phase)));
+ realPhaseOptFourier->SetBinError(i+1, 0.0);
+ }
+
+ return realPhaseOptFourier;
+}
+
+//--------------------------------------------------------------------------
+// GetImaginaryFourier (public)
//--------------------------------------------------------------------------
/**
* returns the imaginary part Fourier as a histogram.
@@ -301,7 +402,7 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
-// GetPowerFourier
+// GetPowerFourier (public)
//--------------------------------------------------------------------------
/**
*
returns the Fourier power spectrum as a histogram.
@@ -343,7 +444,7 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
-// GetPhaseFourier
+// GetPhaseFourier (public)
//--------------------------------------------------------------------------
/**
*
returns the Fourier phase spectrum as a histogram.
@@ -402,7 +503,7 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
-// PrepareFFTwInputData
+// PrepareFFTwInputData (private)
//--------------------------------------------------------------------------
/**
*
Feeds the Fourier data and apply the apodization.
@@ -451,7 +552,7 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
}
//--------------------------------------------------------------------------
-// ApodizeData
+// ApodizeData (private)
//--------------------------------------------------------------------------
/**
*
Carries out the appodization of the data.
diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp
index e163fe68..f83da656 100644
--- a/src/classes/PFourierCanvas.cpp
+++ b/src/classes/PFourierCanvas.cpp
@@ -227,6 +227,7 @@ PFourierCanvas::~PFourierCanvas()
delete fFourierHistos[i].dataFourierIm;
delete fFourierHistos[i].dataFourierPwr;
delete fFourierHistos[i].dataFourierPhase;
+ delete fFourierHistos[i].dataFourierPhaseOptReal;
}
fFourierHistos.clear();
}
@@ -386,6 +387,11 @@ void PFourierCanvas::HandleMenuPopup(Int_t id)
fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE);
fCurrentPlotView = FOURIER_PLOT_PHASE;
PlotFourier();
+ } 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();
} 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) {
@@ -590,6 +596,12 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName)
yAxis = TString("");
xMinBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetFirst();
xMaxBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetLast();
+ break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ xAxis = fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetTitle();
+ yAxis = TString("");
+ xMinBin = fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetFirst();
+ xMaxBin = fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetLast();
break;
default:
xAxis = TString("??");
@@ -624,6 +636,12 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName)
xMinBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetFirst();
xMaxBin = fFourierHistos[0].dataFourierPhase->GetXaxis()->GetLast();
break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ xAxis = fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetTitle();
+ yAxis = TString("Phase");
+ xMinBin = fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetFirst();
+ xMaxBin = fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetLast();
+ break;
default:
xAxis = TString("??");
yAxis = TString("??");
@@ -659,6 +677,9 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName)
case FOURIER_PLOT_PHASE:
fout << fFourierAverage[0].dataFourierPhase->GetBinCenter(i) << ", " << fFourierAverage[0].dataFourierPhase->GetBinContent(i) << endl;
break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ fout << fFourierAverage[0].dataFourierPhaseOptReal->GetBinCenter(i) << ", " << fFourierAverage[0].dataFourierPhaseOptReal->GetBinContent(i) << endl;
+ break;
default:
break;
}
@@ -693,6 +714,9 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName)
case FOURIER_PLOT_PHASE:
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;
default:
break;
}
@@ -710,6 +734,9 @@ void PFourierCanvas::ExportData(const Char_t *pathFileName)
case FOURIER_PLOT_PHASE:
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;
default:
break;
}
@@ -777,6 +804,8 @@ 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]);
+ cout << "debug> histo[" << i << "]: opt. phase = " << fFourierHistos[i].optPhase * 180.0 / TMath::Pi() << "°" << endl;
}
// rescale histo to abs(maximum) == 1
@@ -796,6 +825,7 @@ void PFourierCanvas::InitFourierDataSets()
fFourierHistos[i].dataFourierRe->SetBinContent(j, fFourierHistos[i].dataFourierRe->GetBinContent(j)/fabs(max));
}
}
+
// imaginary
for (UInt_t i=0; iSetBinContent(j, fFourierHistos[i].dataFourierIm->GetBinContent(j)/fabs(max));
}
}
+
// power
max = 0.0;
for (UInt_t i=0; iSetBinContent(j, fFourierHistos[i].dataFourierPwr->GetBinContent(j)/fabs(max));
}
}
+
// phase
max = 0.0;
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));
+ }
+ }
+
// set the marker and line color
for (UInt_t i=0; iSetMarkerColor(fColorList[i]);
@@ -848,6 +894,8 @@ 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]);
}
// set the marker symbol and size
@@ -860,6 +908,8 @@ 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);
}
// initialize average histos
@@ -913,6 +963,7 @@ void PFourierCanvas::InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t w
fPopupFourier->AddEntry("Show Real+Imag", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG);
fPopupFourier->AddEntry("Show Power", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR);
fPopupFourier->AddEntry("Show Phase", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE);
+ fPopupFourier->AddEntry("Show Phase Opt Fourier", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_OPT_REAL);
fPopupFourier->AddSeparator();
fPopupFourier->AddEntry("Phase +", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_PLUS);
fPopupFourier->AddEntry("Phase -", P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_MINUS);
@@ -935,6 +986,9 @@ void PFourierCanvas::InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t w
case FOURIER_PLOT_PHASE:
fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE);
break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_OPT_REAL);
+ break;
default:
break;
}
@@ -1029,6 +1083,10 @@ void PFourierCanvas::CleanupAverage()
delete fFourierAverage[i].dataFourierPhase;
fFourierAverage[i].dataFourierPhase = 0;
}
+ if (fFourierAverage[i].dataFourierPhaseOptReal) {
+ delete fFourierAverage[i].dataFourierPhaseOptReal;
+ fFourierAverage[i].dataFourierPhaseOptReal = 0;
+ }
}
fFourierAverage.clear();
}
@@ -1074,6 +1132,11 @@ 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());
+
// real average
for (Int_t j=0; jGetNbinsX(); j++) {
dval = 0.0;
@@ -1133,6 +1196,21 @@ void PFourierCanvas::HandleAverage()
fFourierAverage[0].dataFourierPhase->SetLineColor(kBlack);
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));
+ }
+ 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);
}
// check if per data set shall be averaged
@@ -1196,6 +1274,11 @@ 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());
+
// real average
for (Int_t j=0; jGetNbinsX(); j++) {
dval = 0.0;
@@ -1255,6 +1338,21 @@ void PFourierCanvas::HandleAverage()
fFourierAverage[i].dataFourierPhase->SetLineColor(fColorList[i]);
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));
+ }
+ 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
}
}
}
@@ -1381,6 +1479,26 @@ void PFourierCanvas::PlotFourier()
fFourierHistos[i].dataFourierPhase->Draw("psame");
}
break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax);
+ ymin = GetMinimum(fFourierHistos[0].dataFourierPhaseOptReal, xmin, xmax);
+ ymax = GetMaximum(fFourierHistos[0].dataFourierPhaseOptReal, xmin, xmax);
+ for (UInt_t i=1; i ymax)
+ ymax = GetMaximum(fFourierHistos[i].dataFourierPhaseOptReal, xmin, xmax);
+ if (GetMinimum(fFourierHistos[i].dataFourierPhaseOptReal, xmin, xmax) < ymin)
+ ymin = GetMinimum(fFourierHistos[i].dataFourierPhaseOptReal, xmin, xmax);
+ }
+ fFourierHistos[0].dataFourierPhaseOptReal->GetYaxis()->SetRangeUser(ymin, 1.05*ymax);
+ fFourierHistos[0].dataFourierPhaseOptReal->GetYaxis()->SetTitleOffset(1.3);
+ fFourierHistos[0].dataFourierPhaseOptReal->GetYaxis()->SetDecimals(kTRUE);
+ fFourierHistos[0].dataFourierPhaseOptReal->GetYaxis()->SetTitle("Phase Opt. Real");
+ fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->SetTitle(fXaxisTitle.Data());
+ fFourierHistos[0].dataFourierPhaseOptReal->Draw("p");
+ for (UInt_t i=1; iDraw("psame");
+ }
+ break;
default:
break;
}
@@ -1541,6 +1659,24 @@ void PFourierCanvas::PlotAverage()
fFourierAverage[0].dataFourierPhase->GetXaxis()->SetRangeUser(xmin, xmax);
fFourierAverage[0].dataFourierPhase->Draw("p");
break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ fFourierAverage[0].dataFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax);
+ ymin = GetMinimum(fFourierAverage[0].dataFourierPhaseOptReal, xmin, xmax);
+ ymax = GetMaximum(fFourierAverage[0].dataFourierPhaseOptReal, xmin, xmax);
+ for (UInt_t i=1; i ymax)
+ ymax = GetMaximum(fFourierAverage[i].dataFourierPhaseOptReal, xmin, xmax);
+ if (GetMinimum(fFourierAverage[i].dataFourierPhaseOptReal, xmin, xmax) < ymin)
+ ymin = GetMinimum(fFourierAverage[i].dataFourierPhaseOptReal, xmin, xmax);
+ }
+ fFourierAverage[0].dataFourierPhaseOptReal->GetYaxis()->SetRangeUser(ymin, 1.03*ymax);
+ fFourierAverage[0].dataFourierPhaseOptReal->GetYaxis()->SetTitleOffset(1.3);
+ fFourierAverage[0].dataFourierPhaseOptReal->GetYaxis()->SetDecimals(kTRUE);
+ fFourierAverage[0].dataFourierPhaseOptReal->GetYaxis()->SetTitle("");
+ fFourierAverage[0].dataFourierPhaseOptReal->GetXaxis()->SetTitle(fXaxisTitle.Data());
+ fFourierAverage[0].dataFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax);
+ fFourierAverage[0].dataFourierPhaseOptReal->Draw("p");
+ break;
default:
break;
}
@@ -1568,6 +1704,9 @@ void PFourierCanvas::PlotAverage()
case FOURIER_PLOT_PHASE:
fFourierAverage[i].dataFourierPhase->Draw("psame");
break;
+ case FOURIER_PLOT_PHASE_OPT_REAL:
+ fFourierAverage[i].dataFourierPhaseOptReal->Draw("psame");
+ break;
default:
break;
}
diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp
index 67c59201..8ce87844 100644
--- a/src/classes/PMsrHandler.cpp
+++ b/src/classes/PMsrHandler.cpp
@@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
- * Copyright (C) 2007-2014 by Andreas Suter *
+ * Copyright (C) 2007-2015 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -242,7 +242,7 @@ Int_t PMsrHandler::ReadMsrFile()
if (result == PMUSR_SUCCESS)
if (!HandleFunctionsEntry(functions))
result = PMUSR_MSR_SYNTAX_ERROR;
- if (result == PMUSR_SUCCESS)
+ if ((result == PMUSR_SUCCESS) && (global.size()>0))
if (!HandleGlobalEntry(global))
result = PMUSR_MSR_SYNTAX_ERROR;
if (result == PMUSR_SUCCESS)
@@ -1595,7 +1595,106 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co
// write GLOBAL block
if (fGlobal.IsPresent()) {
- // not sure that anything needs to be done here ...
+ fout << "GLOBAL" << endl;
+
+ // fittype
+ if (fGlobal.GetFitType() != -1) {
+ fout.width(16);
+ switch (fGlobal.GetFitType()) {
+ case MSR_FITTYPE_SINGLE_HISTO:
+ fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
+ break;
+ case MSR_FITTYPE_ASYM:
+ fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
+ break;
+ case MSR_FITTYPE_MU_MINUS:
+ fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
+ break;
+ case MSR_FITTYPE_NON_MUSR:
+ fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ;
+ break;
+ default:
+ break;
+ }
+ }
+
+ // data range
+ if ((fGlobal.GetDataRange(0) != -1) || (fGlobal.GetDataRange(1) != -1) || (fGlobal.GetDataRange(2) != -1) || (fGlobal.GetDataRange(3) != -1)) {
+ fout.width(16);
+ fout << left << "data";
+ for (UInt_t j=0; j<4; ++j) {
+ if (fGlobal.GetDataRange(j) > 0) {
+ fout.width(8);
+ fout << left << fGlobal.GetDataRange(j);
+ }
+ }
+ fout << endl;
+ }
+
+ // t0
+ if (fGlobal.GetT0BinSize() > 0) {
+ fout.width(16);
+ fout << left << "t0";
+ for (UInt_t j=0; j 0) {
+ fout.width(16);
+ fout << left << "addt0";
+ for (Int_t k=0; k 0)
+ fout << "+" << fGlobal.GetFitRangeOffset(0);
+ fout << " lgb";
+ if (fGlobal.GetFitRangeOffset(1) > 0)
+ fout << "-" << fGlobal.GetFitRangeOffset(1);
+ } else { // fit range given in time
+ for (UInt_t j=0; j<2; j++) {
+ if (fGlobal.GetFitRange(j) == -1)
+ break;
+ UInt_t neededWidth = 7;
+ UInt_t neededPrec = LastSignificant(fRuns[i].GetFitRange(j));
+ fout.width(neededWidth);
+ fout.precision(neededPrec);
+ fout << left << fixed << fGlobal.GetFitRange(j);
+ if (j==0)
+ fout << " ";
+ }
+ }
+ fout << endl;
+ }
+
+ // packing
+ if (fGlobal.GetPacking() != -1) {
+ fout.width(16);
+ fout << left << "packing";
+ fout << fGlobal.GetPacking() << endl;
+ }
+
+ fout << endl << hline.Data() << endl;
}
// write RUN blocks
@@ -1660,22 +1759,24 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co
}
// fittype
- fout.width(16);
- switch (fRuns[i].GetFitType()) {
- case MSR_FITTYPE_SINGLE_HISTO:
- fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
- break;
- case MSR_FITTYPE_ASYM:
- fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
- break;
- case MSR_FITTYPE_MU_MINUS:
- fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
- break;
- case MSR_FITTYPE_NON_MUSR:
- fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ;
- break;
- default:
- break;
+ if (fRuns[i].GetFitType() != -1) {
+ fout.width(16);
+ switch (fRuns[i].GetFitType()) {
+ case MSR_FITTYPE_SINGLE_HISTO:
+ fout << left << "fittype" << MSR_FITTYPE_SINGLE_HISTO << " (single histogram fit)" << endl;
+ break;
+ case MSR_FITTYPE_ASYM:
+ fout << left << "fittype" << MSR_FITTYPE_ASYM << " (asymmetry fit)" << endl ;
+ break;
+ case MSR_FITTYPE_MU_MINUS:
+ fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ;
+ break;
+ case MSR_FITTYPE_NON_MUSR:
+ fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ;
+ break;
+ default:
+ break;
+ }
}
// alpha
@@ -1815,17 +1916,19 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co
}
// addt0
- for (UInt_t j = 0; j < fRuns[i].GetRunNameSize() - 1; ++j) {
- if (fRuns[i].GetAddT0BinSize(j) > 0) {
- fout.width(16);
- fout << left << "addt0";
- for (Int_t k=0; k 0) {
+ for (UInt_t j = 0; j < fRuns[i].GetRunNameSize() - 1; ++j) {
+ if (fRuns[i].GetAddT0BinSize(j) > 0) {
+ fout.width(16);
+ fout << left << "addt0";
+ for (Int_t k=0; k *co
}
// fit
- fout.width(16);
- fout << left << "fit";
- if (fRuns[i].IsFitRangeInBin()) { // fit range given in bins
- fout << "fgb";
- if (fRuns[i].GetFitRangeOffset(0) > 0)
- fout << "+" << fRuns[i].GetFitRangeOffset(0);
- fout << " lgb";
- if (fRuns[i].GetFitRangeOffset(1) > 0)
- fout << "-" << fRuns[i].GetFitRangeOffset(1);
- } else { // fit range given in time
- for (UInt_t j=0; j<2; j++) {
- if (fRuns[i].GetFitRange(j) == -1)
- break;
- UInt_t neededWidth = 7;
- UInt_t neededPrec = LastSignificant(fRuns[i].GetFitRange(j));
- fout.width(neededWidth);
- fout.precision(neededPrec);
- fout << left << fixed << fRuns[i].GetFitRange(j);
- if (j==0)
- fout << " ";
+ if ( (fRuns[i].IsFitRangeInBin() && fRuns[i].GetFitRangeOffset(0) != -1) ||
+ (fRuns[i].GetFitRange(0) != PMUSR_UNDEFINED) ) {
+ fout.width(16);
+ fout << left << "fit";
+ if (fRuns[i].IsFitRangeInBin()) { // fit range given in bins
+ fout << "fgb";
+ if (fRuns[i].GetFitRangeOffset(0) > 0)
+ fout << "+" << fRuns[i].GetFitRangeOffset(0);
+ fout << " lgb";
+ if (fRuns[i].GetFitRangeOffset(1) > 0)
+ fout << "-" << fRuns[i].GetFitRangeOffset(1);
+ } else { // fit range given in time
+ for (UInt_t j=0; j<2; j++) {
+ if (fRuns[i].GetFitRange(j) == -1)
+ break;
+ UInt_t neededWidth = 7;
+ UInt_t neededPrec = LastSignificant(fRuns[i].GetFitRange(j));
+ fout.width(neededWidth);
+ fout.precision(neededPrec);
+ fout << left << fixed << fRuns[i].GetFitRange(j);
+ if (j==0)
+ fout << " ";
+ }
}
+ fout << endl;
}
- fout << endl;
// packing
- fout.width(16);
- fout << left << "packing";
- fout << fRuns[i].GetPacking() << endl;
+ if (fRuns[i].GetPacking() != -1) {
+ fout.width(16);
+ fout << left << "packing";
+ fout << fRuns[i].GetPacking() << endl;
+ }
+
fout << endl;
}
diff --git a/src/include/PFourier.h b/src/include/PFourier.h
index 868b1d09..3e120b7c 100644
--- a/src/include/PFourier.h
+++ b/src/include/PFourier.h
@@ -57,6 +57,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* 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);
diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h
index 634abde4..096097c1 100644
--- a/src/include/PFourierCanvas.h
+++ b/src/include/PFourierCanvas.h
@@ -43,23 +43,26 @@
#define P_MENU_ID_AVERAGE_PER_DATA_SET 10003
#define P_MENU_ID_EXPORT_DATA 10004
-#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
//------------------------------------------------------------------------
/**
* Structure holding all necessary Fourier histograms.
*/
typedef struct {
- TH1F *dataFourierRe; ///< real part of the Fourier transform of the data histogram
- 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 *dataFourierRe; ///< real part of the Fourier transform of the data histogram
+ 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 otpimized real Fourier transform of the data histogram
+ Double_t optPhase; ///< optimal phase which maximizes the real Fourier
} PFourierCanvasDataSet;
//------------------------------------------------------------------------
diff --git a/src/include/PMusr.h b/src/include/PMusr.h
index 8d2683c3..9074a34b 100644
--- a/src/include/PMusr.h
+++ b/src/include/PMusr.h
@@ -123,12 +123,13 @@ typedef struct { char a[7]; } __float128; // needed since cint doesn't know it
#define FOURIER_APOD_MEDIUM 3
#define FOURIER_APOD_STRONG 4
-#define FOURIER_PLOT_NOT_GIVEN 0
-#define FOURIER_PLOT_REAL 1
-#define FOURIER_PLOT_IMAG 2
-#define FOURIER_PLOT_REAL_AND_IMAG 3
-#define FOURIER_PLOT_POWER 4
-#define FOURIER_PLOT_PHASE 5
+#define FOURIER_PLOT_NOT_GIVEN 0
+#define FOURIER_PLOT_REAL 1
+#define FOURIER_PLOT_IMAG 2
+#define FOURIER_PLOT_REAL_AND_IMAG 3
+#define FOURIER_PLOT_POWER 4
+#define FOURIER_PLOT_PHASE 5
+#define FOURIER_PLOT_PHASE_OPT_REAL 6
//-------------------------------------------------------------
// RRF related tags
diff --git a/src/musrFT.cpp b/src/musrFT.cpp
index 37c931c3..84bc5991 100644
--- a/src/musrFT.cpp
+++ b/src/musrFT.cpp
@@ -67,7 +67,7 @@ typedef struct {
TString msrFlnOut; ///< dump file name for msr-file generation
Int_t bkg_range[2]; ///< background range
PDoubleVector bkg; ///< background value
- TString fourierOpt; ///< Fourier options, i.e. real, imag, power, phase
+ TString fourierOpt; ///< Fourier options, i.e. real, imag, power, phase, phaseOptReal
TString apodization; ///< apodization setting: none, weak, medium, strong
Int_t fourierPower; ///< Fourier power for zero padding, i.e. 2^fourierPower points
TString fourierUnits; ///< wished Fourier units: Gauss, Tesla, MHz, Mc/s
@@ -111,7 +111,7 @@ void musrFT_syntax()
cout << endl << " -br, --background-range : background interval used to estimate the background to be";
cout << endl << " subtracted before the Fourier transform. , to be given in bins.";
cout << endl << " -bg, --background : gives the background explicit for each histogram.";
- cout << endl << " -fo, --fourier-option : can be 'real', 'imag', 'real+imag', 'power', or 'phase'.";
+ cout << endl << " -fo, --fourier-option : can be 'real', 'imag', 'real+imag', 'power', 'phase', or 'phaseOptReal'.";
cout << endl << " If this is not defined (neither on the command line nor in the musrfit_startup.xml),";
cout << endl << " default will be 'power'.";
cout << endl << " -ap, --apodization : can be either 'none', 'weak', 'medium', 'strong'.";
@@ -370,7 +370,8 @@ Int_t musrFT_parse_options(Int_t argc, Char_t *argv[], musrFT_startup_param &sta
return 2;
}
TString topt(argv[i+1]);
- if (!topt.BeginsWith("real") && !topt.BeginsWith("imag") && !topt.BeginsWith("power") && !topt.BeginsWith("phase")) {
+ if (!topt.BeginsWith("real") && !topt.BeginsWith("imag") && !topt.BeginsWith("power") &&
+ !topt.BeginsWith("phase") && !topt.BeginsWith("phaseOptReal")) {
cerr << endl << ">> musrFT **ERROR** found option --fourier-option with unrecognized argument '" << topt << "'." << endl;
return 2;
}
@@ -941,9 +942,9 @@ void musrFT_dumpMsrFile(musrFT_startup_param ¶m)
fout << "units " << param.fourierUnits << " # units either 'Gauss', 'MHz', or 'Mc/s'" << endl;
}
if (param.fourierOpt.BeginsWith("??")) { // Fourier plot option not given, hence choose POWER
- fout << "plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE" << endl;
+ fout << "plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL" << endl;
} else {
- fout << "plot " << param.fourierOpt << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE" << endl;
+ fout << "plot " << param.fourierOpt << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL" << endl;
}
if (param.fourierPower > 1) {
fout << "fourier_power " << param.fourierPower << endl;
@@ -1368,6 +1369,8 @@ Int_t main(Int_t argc, Char_t *argv[])
fourierPlotTag = FOURIER_PLOT_POWER;
else if (!startupParam.fourierOpt.CompareTo("phase", TString::kIgnoreCase))
fourierPlotTag = FOURIER_PLOT_PHASE;
+ else if (!startupParam.fourierOpt.CompareTo("phaseoptreal", TString::kIgnoreCase))
+ fourierPlotTag = FOURIER_PLOT_PHASE_OPT_REAL;
else
fourierPlotTag = FOURIER_PLOT_POWER;
}
diff --git a/src/musredit/PGetMusrFTOptionsDialog.cpp b/src/musredit/PGetMusrFTOptionsDialog.cpp
index 5be9bc08..cdf15e3f 100644
--- a/src/musredit/PGetMusrFTOptionsDialog.cpp
+++ b/src/musredit/PGetMusrFTOptionsDialog.cpp
@@ -36,12 +36,13 @@
#include "PGetMusrFTOptionsDialog.h"
-#define MUSRFT_OPT_UNDEF 0
-#define MUSRFT_OPT_REAL 1
-#define MUSRFT_OPT_IMAG 2
-#define MUSRFT_OPT_REAL_AND_IMAG 3
-#define MUSRFT_OPT_POWER 4
-#define MUSRFT_OPT_PHASE 5
+#define MUSRFT_OPT_UNDEF 0
+#define MUSRFT_OPT_REAL 1
+#define MUSRFT_OPT_IMAG 2
+#define MUSRFT_OPT_REAL_AND_IMAG 3
+#define MUSRFT_OPT_POWER 4
+#define MUSRFT_OPT_PHASE 5
+#define MUSRFT_OPT_PHASE_OPT_REAL 6
#define MUSRFT_APOD_UNDEF 0
#define MUSRFT_APOD_WEAK 1
@@ -139,6 +140,8 @@ PGetMusrFTOptionsDialog::PGetMusrFTOptionsDialog(QString currentMsrFile, QString
fFourierOption_comboBox->setCurrentIndex(MUSRFT_OPT_POWER);
else if (prevCmd[i+1] == "phase")
fFourierOption_comboBox->setCurrentIndex(MUSRFT_OPT_PHASE);
+ else if (prevCmd[i+1] == "phaseOptReal")
+ fFourierOption_comboBox->setCurrentIndex(MUSRFT_OPT_PHASE_OPT_REAL);
else
fFourierOption_comboBox->setCurrentIndex(MUSRFT_OPT_UNDEF);
i++;
diff --git a/src/musredit/forms/PGetMusrFTOptionsDialog.ui b/src/musredit/forms/PGetMusrFTOptionsDialog.ui
index 75966be9..6935f643 100644
--- a/src/musredit/forms/PGetMusrFTOptionsDialog.ui
+++ b/src/musredit/forms/PGetMusrFTOptionsDialog.ui
@@ -34,7 +34,7 @@
Fourier
-
+
21
@@ -91,6 +91,11 @@ p, li { white-space: pre-wrap; }
phase
+ -
+
+ phaseOptReal
+
+