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

first almost feature complete version of musrFT

Browse Source
This commit is contained in:
suter_a
2015-02-13 16:56:19 +01:00
parent d4ba2f6d81
commit 2d2594c903
12 changed files with 1900 additions and 275 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/classes/Makefile.am
View File

@@ -4,6 +4,7 @@ h_sources = \
../include/PFitterFcn.h \
../include/PFitter.h \
../include/PFourier.h \
../include/PFourierCanvas.h \
../include/PFunctionGrammar.h \
../include/PFunction.h \
../include/PFunctionHandler.h \
@@ -27,6 +28,7 @@ h_sources_userFcn = \
../include/PUserFcnBase.h
h_linkdef = \
../include/PFourierCanvasLinkDef.h \
../include/PMusrCanvasLinkDef.h \
../include/PMusrT0LinkDef.h \
../include/PStartupHandlerLinkDef.h
@@ -35,6 +37,7 @@ h_linkdef_userFcn = \
../include/PUserFcnBaseLinkDef.h
dict_h_sources = \
PFourierCanvasDict.h \
PMusrCanvasDict.h \
PMusrT0Dict.h \
PStartupHandlerDict.h
@@ -46,6 +49,7 @@ cpp_sources = \
PFitter.cpp \
PFitterFcn.cpp \
PFourier.cpp \
PFourierCanvas.cpp \
PFunction.cpp \
PFunctionHandler.cpp \
PMsr2Data.cpp \
@@ -68,6 +72,7 @@ cpp_sources_userFcn = \
PUserFcnBase.cpp
dict_cpp_sources = \
PFourierCanvasDict.cpp \
PMusrCanvasDict.cpp \
PMusrT0Dict.cpp \
PStartupHandlerDict.cpp

1312
src/classes/PFourierCanvas.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

203
src/classes/PMusrCanvas.cpp
View File

@@ -3797,15 +3797,10 @@ void PMusrCanvas::HandleAverage()
// calculate all the average data sets
double dval;
if (fDataAvg.data != 0) {
if (!CalcAlignment(eTime)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: data: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].data->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].data->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].data, fData[0].data->GetBinContent(i));
}
fDataAvg.data->SetBinContent(i, dval/fData.size());
}
@@ -3816,15 +3811,10 @@ void PMusrCanvas::HandleAverage()
fDataAvg.data->SetMarkerStyle(fData[0].data->GetMarkerStyle());
}
if (fDataAvg.dataFourierRe != 0) {
if (!CalcAlignment(eFreq)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Re: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].dataFourierRe->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].dataFourierRe->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].dataFourierRe, fData[0].dataFourierRe->GetBinContent(i));
}
fDataAvg.dataFourierRe->SetBinContent(i, dval/fData.size());
}
@@ -3835,15 +3825,10 @@ void PMusrCanvas::HandleAverage()
fDataAvg.dataFourierRe->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle());
}
if (fDataAvg.dataFourierIm != 0) {
if (!CalcAlignment(eFreq)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Im: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].dataFourierIm->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].dataFourierIm->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].dataFourierIm, fData[0].dataFourierIm->GetBinContent(i));
}
fDataAvg.dataFourierIm->SetBinContent(i, dval/fData.size());
}
@@ -3854,15 +3839,10 @@ void PMusrCanvas::HandleAverage()
fDataAvg.dataFourierIm->SetMarkerStyle(fData[0].dataFourierIm->GetMarkerStyle());
}
if (fDataAvg.dataFourierPwr != 0) {
if (!CalcAlignment(eFreq)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Pwr: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].dataFourierPwr->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].dataFourierPwr->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].dataFourierPwr, fData[0].dataFourierPwr->GetBinContent(i));
}
fDataAvg.dataFourierPwr->SetBinContent(i, dval/fData.size());
}
@@ -3873,15 +3853,10 @@ void PMusrCanvas::HandleAverage()
fDataAvg.dataFourierPwr->SetMarkerStyle(fData[0].dataFourierPwr->GetMarkerStyle());
}
if (fDataAvg.dataFourierPhase != 0) {
if (!CalcAlignment(eFreq)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Phase: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].dataFourierPhase->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].dataFourierPhase->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].dataFourierPhase, fData[0].dataFourierPhase->GetBinContent(i));
}
fDataAvg.dataFourierPhase->SetBinContent(i, dval/fData.size());
}
@@ -3892,15 +3867,10 @@ void PMusrCanvas::HandleAverage()
fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle());
}
if (fDataAvg.theory != 0) {
if (!CalcAlignment(eTheoTime)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: theory: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].theory->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].theory->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].theory, fData[0].theory->GetBinContent(i));
}
fDataAvg.theory->SetBinContent(i, dval/fData.size());
}
@@ -3910,7 +3880,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].theoryFourierRe->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].theoryFourierRe->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].theoryFourierRe, fData[0].theoryFourierRe->GetBinContent(i));
}
fDataAvg.theoryFourierRe->SetBinContent(i, dval/fData.size());
}
@@ -3924,7 +3894,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].theoryFourierIm->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].theoryFourierIm->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].theoryFourierIm, fData[0].theoryFourierIm->GetBinContent(i));
}
fDataAvg.theoryFourierIm->SetBinContent(i, dval/fData.size());
}
@@ -3938,7 +3908,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].theoryFourierPwr->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].theoryFourierPwr->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].theoryFourierPwr, fData[0].theoryFourierPwr->GetBinContent(i));
}
fDataAvg.theoryFourierPwr->SetBinContent(i, dval/fData.size());
}
@@ -3952,7 +3922,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].theoryFourierPhase->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].theoryFourierPhase->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].theoryFourierPhase, fData[0].theoryFourierPhase->GetBinContent(i));
}
fDataAvg.theoryFourierPhase->SetBinContent(i, dval/fData.size());
}
@@ -3963,15 +3933,10 @@ void PMusrCanvas::HandleAverage()
fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle());
}
if (fDataAvg.diff != 0) {
if (!CalcAlignment(eTime)) {
cerr << endl << ">> PMusrCanvas::HandleAverage: diff: **WARNING** only approx. alignment possible." << endl;
}
for (Int_t i=0; i<fData[0].diff->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
if ((i-fAlignmentOffset[j]) > 0) {
dval += fData[j].diff->GetBinContent(i-fAlignmentOffset[j]);
}
dval += GetInterpolatedValue(fData[j].diff, fData[0].diff->GetBinContent(i));
}
fDataAvg.diff->SetBinContent(i, dval/fData.size());
}
@@ -3985,7 +3950,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].diffFourierRe->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].diffFourierRe->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].diffFourierRe, fData[0].diffFourierRe->GetBinContent(i));
}
fDataAvg.diffFourierRe->SetBinContent(i, dval/fData.size());
}
@@ -3999,7 +3964,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].diffFourierIm->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].diffFourierIm->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].diffFourierIm, fData[0].diffFourierIm->GetBinContent(i));
}
fDataAvg.diffFourierIm->SetBinContent(i, dval/fData.size());
}
@@ -4013,7 +3978,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].diffFourierPwr->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].diffFourierPwr->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].diffFourierPwr, fData[0].diffFourierPwr->GetBinContent(i));
}
fDataAvg.diffFourierPwr->SetBinContent(i, dval/fData.size());
}
@@ -4027,7 +3992,7 @@ void PMusrCanvas::HandleAverage()
for (Int_t i=0; i<fData[0].diffFourierPhase->GetNbinsX(); i++) {
dval = 0.0;
for (UInt_t j=0; j<fData.size(); j++) {
dval += fData[j].diffFourierPhase->GetBinContent(i);
dval += GetInterpolatedValue(fData[j].diffFourierPhase, fData[0].diffFourierPhase->GetBinContent(i));
}
fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size());
}
@@ -6314,130 +6279,40 @@ UInt_t PMusrCanvas::GetNeededAccuracy(PMsrParamStructure param)
//--------------------------------------------------------------------------
// CalcAlignment (private)
// GetInterpolatedValue (private)
//--------------------------------------------------------------------------
/**
* <p>Calculates the alignment index for each data set needed to average the data.
* <p>search for xVal in histo. If xVal is not found exactly, interpolate and
* return the interpolated y-value.
*
* <b>return:</b>
* - true for perfect alignment
* - false for approximate alignment
* - interpolated value if xVal is within histo range, 0 otherwise.
*
* \param tag to distinguish time data sets from Fourier data sets.
* \param histo pointer of the histogram
* \param xVal x-value to be looked for
*/
Bool_t PMusrCanvas::CalcAlignment(const EAlignTag tag)
Double_t PMusrCanvas::GetInterpolatedValue(TH1F* histo, Double_t xVal)
{
Bool_t result = true;
if (histo == 0)
return 0.0;
fAlignmentOffset.clear();
Int_t idx = histo->FindBin(xVal);
UInt_t idx=0;
Double_t dval;
if (tag == eTime) {
// first find the data vector with the lowest initial time
dval = fData[0].data->GetXaxis()->GetBinCenter(1);
for (UInt_t i=1; i<fData.size(); i++) {
if (fData[i].data->GetXaxis()->GetBinCenter(1) < dval) {
idx = i;
dval = fData[i].data->GetXaxis()->GetBinCenter(1);
}
}
// make sure idx is within range
if ((idx < 1) || (idx > histo->GetNbinsX()))
return 0.0;
// next setp: find all the alignment indices
fAlignmentOffset.resize(fData.size());
for (UInt_t i=0; i<fData.size(); i++) {
dval = fData[i].data->GetXaxis()->GetBinCenter(1);
Int_t j=1;
Bool_t found = false;
do {
if (fData[idx].data->GetXaxis()->GetBinCenter(j) >= dval) {
fAlignmentOffset[i] = (UInt_t)(j-1);
if (fData[idx].data->GetXaxis()->GetBinCenter(j) != dval) {
result = false;
}
found = true;
}
} while (!found && (++j<fData[idx].data->GetNbinsX()));
}
} else if (tag == eTheoTime) {
// first find the data vector with the lowest initial time
dval = fData[0].theory->GetXaxis()->GetBinCenter(1);
for (UInt_t i=1; i<fData.size(); i++) {
if (fData[i].theory->GetXaxis()->GetBinCenter(1) < dval) {
idx = i;
dval = fData[i].theory->GetXaxis()->GetBinCenter(1);
}
}
// make sure the lower bound idx is found. This is needed since
// FindBin rounds towards the closer idx.
if (histo->GetBinCenter(idx) > xVal)
--idx;
// next setp: find all the alignment indices
fAlignmentOffset.resize(fData.size());
for (UInt_t i=0; i<fData.size(); i++) {
dval = fData[i].theory->GetXaxis()->GetBinCenter(1);
Int_t j=1;
Bool_t found = false;
do {
if (fData[idx].theory->GetXaxis()->GetBinCenter(j) >= dval) {
fAlignmentOffset[i] = (UInt_t)(j-1);
if (fData[idx].theory->GetXaxis()->GetBinCenter(j) != dval) {
result = false;
}
found = true;
}
} while (!found && (++j<fData[idx].theory->GetNbinsX()));
}
} else if (tag == eFreq) {
// first find the data vector with the lowest initial time
dval = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(1);
for (UInt_t i=1; i<fData.size(); i++) {
if (fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1) < dval) {
idx = i;
dval = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1);
}
}
Double_t x0, x1, y0, y1;
x0 = histo->GetBinCenter(idx);
x1 = histo->GetBinCenter(idx+1);
y0 = histo->GetBinContent(idx);
y1 = histo->GetBinContent(idx+1);
// next setp: find all the alignment indices
fAlignmentOffset.resize(fData.size());
for (UInt_t i=0; i<fData.size(); i++) {
dval = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1);
Int_t j=1;
Bool_t found = false;
do {
if (fData[idx].dataFourierRe->GetXaxis()->GetBinCenter(j) >= dval) {
fAlignmentOffset[i] = (UInt_t)(j-1);
if (fData[idx].dataFourierRe->GetXaxis()->GetBinCenter(j) != dval) {
result = false;
}
found = true;
}
} while (!found && (++j<fData[idx].dataFourierRe->GetNbinsX()));
}
} else if (tag == eTheoFreq) {
// first find the data vector with the lowest initial time
dval = fData[0].theoryFourierRe->GetXaxis()->GetBinCenter(1);
for (UInt_t i=1; i<fData.size(); i++) {
if (fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1) < dval) {
idx = i;
dval = fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1);
}
}
// next setp: find all the alignment indices
fAlignmentOffset.resize(fData.size());
for (UInt_t i=0; i<fData.size(); i++) {
dval = fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1);
Int_t j=1;
Bool_t found = false;
do {
if (fData[idx].theoryFourierRe->GetXaxis()->GetBinCenter(j) >= dval) {
fAlignmentOffset[i] = (UInt_t)(j-1);
if (fData[idx].theoryFourierRe->GetXaxis()->GetBinCenter(j) != dval) {
result = false;
}
found = true;
}
} while (!found && (++j<fData[idx].theoryFourierRe->GetNbinsX()));
}
}
return result;
return (y1-y0)*(xVal-x0)/(x1-x0)+y0;
}

53
src/classes/PPrepFourier.cpp
View File

@@ -70,14 +70,12 @@ PPrepFourier::~PPrepFourier()
// SetBkgRange
//--------------------------------------------------------------------------
/**
* <p>
* <p>set the background range.
*
* \param bkgRange
* \param bkgRange array with background range
*/
void PPrepFourier::SetBkgRange(const Int_t *bkgRange)
{
cout << endl << "debug> bkgRange: " << bkgRange[0] << ", " << bkgRange[1] << endl;
int err=0;
if (bkgRange[0] >= -1) {
fBkgRange[0] = bkgRange[0];
@@ -118,9 +116,9 @@ void PPrepFourier::SetBkgRange(const Int_t *bkgRange)
// SetPacking
//--------------------------------------------------------------------------
/**
* <p>
* <p>set the packing for the histograms.
*
* \param packing
* \param packing number to be used.
*/
void PPrepFourier::SetPacking(const Int_t packing)
{
@@ -135,9 +133,10 @@ void PPrepFourier::SetPacking(const Int_t packing)
// AddData
//--------------------------------------------------------------------------
/**
* <p>
* <p>add a data-set (time domain data + meta information) to the internal
* data vector.
*
* \param data
* \param data set to be added
*/
void PPrepFourier::AddData(musrFT_data &data)
{
@@ -148,13 +147,10 @@ void PPrepFourier::AddData(musrFT_data &data)
// DoBkgCorrection
//--------------------------------------------------------------------------
/**
* <p>
*
* <p>Correct the internal data sets according to a background interval given.
*/
void PPrepFourier::DoBkgCorrection()
{
cout << endl << "debug> DoBkgCorrection ...";
// make sure fData are already present, and if not create the necessary data sets
if (fData.size() != fRawData.size()) {
InitData();
@@ -162,7 +158,6 @@ void PPrepFourier::DoBkgCorrection()
// if no bkg-range is given, nothing needs to be done
if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) {
cout << endl << "debug> no background range given ...";
return;
}
@@ -181,7 +176,6 @@ void PPrepFourier::DoBkgCorrection()
bkg += fRawData[i].rawData[j];
}
bkg /= (fBkgRange[1]-fBkgRange[0]+1);
cout << endl << "debug> histo: " << i << ", bkg=" << bkg;
// correct data
for (unsigned int j=0; j<fData[i].size(); j++)
@@ -193,13 +187,10 @@ void PPrepFourier::DoBkgCorrection()
// DoPacking
//--------------------------------------------------------------------------
/**
* <p>
*
* <p>Rebin (pack) the internal data.
*/
void PPrepFourier::DoPacking()
{
cout << endl << "debug> DoPacking : pack=" << fPacking << " ...";
// make sure fData are already present, and if not create the necessary data sets
if (fData.size() != fRawData.size()) {
InitData();
@@ -224,10 +215,6 @@ void PPrepFourier::DoPacking()
// change the original data set with the packed one
fData[i].clear();
fData[i] = tmpData;
cout << endl << "debug> histo " << i+1 << ": packed data: ";
for (unsigned int j=0; j<15; j++)
cout << fData[i][j] << ", ";
}
}
@@ -235,27 +222,25 @@ void PPrepFourier::DoPacking()
// DoFiltering
//--------------------------------------------------------------------------
/**
* <p>
*
* <p>Not implemented yet.
*/
void PPrepFourier::DoFiltering()
{
cout << endl << "debug> DoFiltering not yet implemented ...";
// make sure fData are already present, and if not create the necessary data sets
if (fData.size() != fRawData.size()) {
InitData();
}
}
//--------------------------------------------------------------------------
// DoLifeTimeCorrection
//--------------------------------------------------------------------------
/**
* <p>
* <p>Try to do a muon life time correction. The idea is to estimate N0 without
* any theory. This will be OK for high fields (> couple kGauss) but not so good
* for low fields.
*
* \param fudge
* \param fudge rescaling factor for the estimated N0. Should be around 1
*/
void PPrepFourier::DoLifeTimeCorrection(Double_t fudge)
{
@@ -294,9 +279,9 @@ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge)
// GetInfo
//--------------------------------------------------------------------------
/**
* <p>
* <p>Returns the meta information of a data set.
*
* \param idx
* \param idx index of the object
*/
TString PPrepFourier::GetInfo(const UInt_t idx)
{
@@ -312,8 +297,8 @@ TString PPrepFourier::GetInfo(const UInt_t idx)
// GetData
//--------------------------------------------------------------------------
/**
* <p>
*
* <p>Creates the requested TH1F objects and returns them. The ownership is with
* the caller.
*/
vector<TH1F*> PPrepFourier::GetData()
{
@@ -382,7 +367,7 @@ vector<TH1F*> PPrepFourier::GetData()
* <p>Creates the requested TH1F object and returns it. The ownership is with
* the caller.
*
* \param idx
* \param idx index of the requested histogram
*/
TH1F *PPrepFourier::GetData(const UInt_t idx)
{