LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

Merge branch 'master' into root6

Browse Source
This commit is contained in:
suter_a 2015-09-24 20:56:04 +02:00
commit 1f7c4898dc
10 changed files with 463 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
ChangeLog
View File

@ -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.

119
src/classes/PFourier.cpp
View File

@ -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)
//--------------------------------------------------------------------------
/**
* <p>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)
//--------------------------------------------------------------------------
/**
* <p>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)
//--------------------------------------------------------------------------
/**
* <p>returns the real part Fourier as a histogram.
@ -259,7 +259,108 @@ TH1F* PFourier::GetRealFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
// GetImaginaryFourier
// GetPhaseOptRealFourier (public)
//--------------------------------------------------------------------------
/**
* <p>returns the phase corrected real Fourier transform. The correction angle is
* past back as well.
* <p>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<Double_t> 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; j<realF.size(); j++) {
sum += realF[j]*cos(da*i) + imagF[j]*sin(da*i);
}
rotRealIntegral[i] = sum;
}
// find the maximum in rotRealIntegral
Double_t maxRot = 0.0;
UInt_t maxRotBin = 0;
for (UInt_t i=0; i<720; i++) {
if (maxRot < rotRealIntegral[i]) {
maxRot = rotRealIntegral[i];
maxRotBin = i;
}
}
// keep the optimal phase
phase = maxRotBin*da;
// clean up
realF.clear();
imagF.clear();
// 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());
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; i<noOfFourierBins; i++) {
realPhaseOptFourier->SetBinContent(i+1, scale*(fOut[i][0]*cos(phase) + fOut[i][1]*sin(phase)));
realPhaseOptFourier->SetBinError(i+1, 0.0);
}
return realPhaseOptFourier;
}
//--------------------------------------------------------------------------
// GetImaginaryFourier (public)
//--------------------------------------------------------------------------
/**
* <p>returns the imaginary part Fourier as a histogram.
@ -301,7 +402,7 @@ TH1F* PFourier::GetImaginaryFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
// GetPowerFourier
// GetPowerFourier (public)
//--------------------------------------------------------------------------
/**
* <p>returns the Fourier power spectrum as a histogram.
@ -343,7 +444,7 @@ TH1F* PFourier::GetPowerFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
// GetPhaseFourier
// GetPhaseFourier (public)
//--------------------------------------------------------------------------
/**
* <p>returns the Fourier phase spectrum as a histogram.
@ -402,7 +503,7 @@ TH1F* PFourier::GetPhaseFourier(const Double_t scale)
}
//--------------------------------------------------------------------------
// PrepareFFTwInputData
// PrepareFFTwInputData (private)
//--------------------------------------------------------------------------
/**
* <p>Feeds the Fourier data and apply the apodization.
@ -451,7 +552,7 @@ void PFourier::PrepareFFTwInputData(UInt_t apodizationTag)
}
//--------------------------------------------------------------------------
// ApodizeData
// ApodizeData (private)
//--------------------------------------------------------------------------
/**
* <p>Carries out the appodization of the data.

139
src/classes/PFourierCanvas.cpp
View File

@ -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) {
@ -591,6 +597,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("??");
@ -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; i<fFourierHistos.size(); i++) {
for (Int_t j=start; j<=end; j++) {
@ -809,6 +839,7 @@ void PFourierCanvas::InitFourierDataSets()
fFourierHistos[i].dataFourierIm->SetBinContent(j, fFourierHistos[i].dataFourierIm->GetBinContent(j)/fabs(max));
}
}
// power
max = 0.0;
for (UInt_t i=0; i<fFourierHistos.size(); i++) {
@ -823,6 +854,7 @@ void PFourierCanvas::InitFourierDataSets()
fFourierHistos[i].dataFourierPwr->SetBinContent(j, fFourierHistos[i].dataFourierPwr->GetBinContent(j)/fabs(max));
}
}
// phase
max = 0.0;
for (UInt_t i=0; i<fFourierHistos.size(); i++) {
@ -838,6 +870,20 @@ void PFourierCanvas::InitFourierDataSets()
}
}
// phase opt real
for (UInt_t i=0; i<fFourierHistos.size(); i++) {
for (Int_t j=start; j<=end; j++) {
dval = fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j);
if (fabs(dval) > max)
max = dval;
}
}
for (UInt_t i=0; i<fFourierHistos.size(); i++) {
for (Int_t j=1; j<fFourierHistos[i].dataFourierPhaseOptReal->GetNbinsX(); j++) {
fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max));
}
}
// set the marker and line color
for (UInt_t i=0; i<fFourierHistos.size(); i++) {
fFourierHistos[i].dataFourierRe->SetMarkerColor(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; j<fFourierHistos[0].dataFourierRe->GetNbinsX(); 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; j<fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(); j++) {
dval = 0.0;
for (UInt_t i=0; i<fFourierHistos.size(); i++) {
if (j < fFourierHistos[i].dataFourierPhaseOptReal->GetNbinsX())
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; j<fFourierHistos[0].dataFourierRe->GetNbinsX(); 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; j<fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(); 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<fFourierHistos.size(); i++) {
if (GetMaximum(fFourierHistos[i].dataFourierPhaseOptReal, xmin, xmax) > 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; i<fFourierHistos.size(); i++) {
fFourierHistos[i].dataFourierPhaseOptReal->Draw("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<fFourierAverage.size(); i++) {
if (GetMaximum(fFourierAverage[i].dataFourierPhaseOptReal, xmin, xmax) > 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("<Phase Opt. Real>");
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;
}

115
src/classes/PMsrHandler.cpp
View File

@ -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<UInt_t, TString> *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<fGlobal.GetT0BinSize(); ++j) {
fout.width(8);
fout.precision(1);
fout.setf(ios::fixed,ios::floatfield);
fout << left << fGlobal.GetT0Bin(j);
}
fout << endl;
}
// addt0
for (UInt_t j = 0; j < fGlobal.GetAddT0BinEntries(); ++j) {
if (fGlobal.GetAddT0BinSize(j) > 0) {
fout.width(16);
fout << left << "addt0";
for (Int_t k=0; k<fGlobal.GetAddT0BinSize(j); ++k) {
fout.width(8);
fout.precision(1);
fout.setf(ios::fixed,ios::floatfield);
fout << left << fGlobal.GetAddT0Bin(j, k);
}
fout << endl;
}
}
// fit range
if ( (fGlobal.IsFitRangeInBin() && fGlobal.GetFitRangeOffset(0) != -1) ||
(fGlobal.GetFitRange(0) != PMUSR_UNDEFINED) ) {
fout.width(16);
fout << left << "fit";
if (fGlobal.IsFitRangeInBin()) { // fit range given in bins
fout << "fgb";
if (fGlobal.GetFitRangeOffset(0) > 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,6 +1759,7 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
}
// fittype
if (fRuns[i].GetFitType() != -1) {
fout.width(16);
switch (fRuns[i].GetFitType()) {
case MSR_FITTYPE_SINGLE_HISTO:
@ -1677,6 +1777,7 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
default:
break;
}
}
// alpha
if (fRuns[i].GetAlphaParamNo() != -1) {
@ -1815,6 +1916,7 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
}
// addt0
if (fRuns[i].GetAddT0BinEntries() > 0) {
for (UInt_t j = 0; j < fRuns[i].GetRunNameSize() - 1; ++j) {
if (fRuns[i].GetAddT0BinSize(j) > 0) {
fout.width(16);
@ -1828,6 +1930,7 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
fout << endl;
}
}
}
// xy-data
if (fRuns[i].GetXDataIndex() != -1) { // indices
@ -1852,6 +1955,8 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
}
// fit
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
@ -1875,11 +1980,15 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
}
}
fout << endl;
}
// packing
if (fRuns[i].GetPacking() != -1) {
fout.width(16);
fout << left << "packing";
fout << fRuns[i].GetPacking() << endl;
}
fout << endl;
}

1
src/include/PFourier.h
View File

@ -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);

7
src/include/PFourierCanvas.h
View File

@ -48,8 +48,9 @@
#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_PHASE_OPT_REAL 105
#define P_MENU_ID_FOURIER_PHASE_PLUS 106
#define P_MENU_ID_FOURIER_PHASE_MINUS 107
//------------------------------------------------------------------------
/**
@ -60,6 +61,8 @@ 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 otpimized real Fourier transform of the data histogram
Double_t optPhase; ///< optimal phase which maximizes the real Fourier
} PFourierCanvasDataSet;
//------------------------------------------------------------------------

1
src/include/PMusr.h
View File

@ -129,6 +129,7 @@ typedef struct { char a[7]; } __float128; // needed since cint doesn't know it
#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

13
src/musrFT.cpp
View File

@ -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 <start> <end>: background interval used to estimate the background to be";
cout << endl << " subtracted before the Fourier transform. <start>, <end> to be given in bins.";
cout << endl << " -bg, --background <list> : gives the background explicit for each histogram.";
cout << endl << " -fo, --fourier-option <fopt>: <fopt> can be 'real', 'imag', 'real+imag', 'power', or 'phase'.";
cout << endl << " -fo, --fourier-option <fopt>: <fopt> 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 <val> : <val> 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 &param)
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;
}

3
src/musredit/PGetMusrFTOptionsDialog.cpp
View File

@ -42,6 +42,7 @@
#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++;

7
src/musredit/forms/PGetMusrFTOptionsDialog.ui
View File

@ -34,7 +34,7 @@
<property name="title">
<string> Fourier </string>
</property>
<widget class="QWidget" name="">
<widget class="QWidget" name="layoutWidget">
<property name="geometry">
<rect>
<x>21</x>
@ -91,6 +91,11 @@ p, li { white-space: pre-wrap; }
<string>phase</string>
</property>
</item>
<item>
<property name="text">
<string>phaseOptReal</string>
</property>
</item>
</widget>
</item>
</layout>