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

added per run block chisq output

Browse Source
This commit is contained in:
nemu 2011-07-25 07:35:29 +00:00
parent 2a4f5bf78a
commit 58a79084b1
18 changed files with 327 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ChangeLog
View File

@ -6,6 +6,7 @@
changes since 0.9.0
===================================
NEW added chisq per run block information (all fit types). The same output rules as for the expected chisq are in place.
NEW calculate expected chisq (Pearson's chisq) for single histogram fits. It always will send this information to the stdout
but only writes it to the msr-file if a corresponding flag in the musrfit_startup.xml is enabling it. The following
expected chisq's are calculated: total expected chisq, and the per histo expected chisq. (MUSR-194)

49
src/classes/PFitter.cpp
View File

@ -174,6 +174,11 @@ Bool_t PFitter::DoFit()
Double_t totalExpectedChisq = 0.0;
std::vector<Double_t> expectedChisqPerHisto;
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
// calculate chisq per run
std::vector<Double_t> chisqPerHisto;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i));
}
cout << endl << endl << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf;
@ -183,14 +188,20 @@ Bool_t PFitter::DoFit()
for (UInt_t i=0; i<expectedChisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i);
if (ndf_histo > 0)
cout << endl << ">> run block " << i+1 << ": expected chisq = " << expectedChisqPerHisto[i] << ", ndf = " << ndf_histo << ", expected chisq/NDF = " << expectedChisqPerHisto[i]/ndf_histo;
else
cout << endl << ">> run block " << i+1 << ": expected chisq = " << expectedChisqPerHisto[i];
cout << endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << "/" << expectedChisqPerHisto[i]/ndf_histo << ")";
}
} else if (chisqPerHisto.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits
UInt_t ndf_histo = 0;
for (UInt_t i=0; i<chisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i);
if (ndf_histo > 0)
cout << endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << ")";
}
// clean up
expectedChisqPerHisto.clear();
}
// clean up
chisqPerHisto.clear();
expectedChisqPerHisto.clear();
} else { // max. log likelihood
cout << endl << endl << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf;
}
@ -1456,6 +1467,12 @@ Bool_t PFitter::ExecuteSave()
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
// calculate chisq per run
std::vector<Double_t> chisqPerHisto;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
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;
@ -1467,11 +1484,31 @@ Bool_t PFitter::ExecuteSave()
// 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; i<chisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(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;

23
src/classes/PFitterFcn.cpp
View File

@ -92,30 +92,11 @@ Double_t PFitterFcn::operator()(const std::vector<Double_t>& par) const
return value;
}
//--------------------------------------------------------------------------
// GetNoOfFittedBins()
//--------------------------------------------------------------------------
/**
* <p>Get the number of fitted bins of the run block idx.
*
* \param idx index of the run block
*/
UInt_t PFitterFcn::GetNoOfFittedBins(const UInt_t idx)
{
UInt_t result = 0;
if (idx < fRunListCollection->GetNoOfSingleHisto())
result = fRunListCollection->GetNoOfBinsFitted(idx);
return result;
}
//--------------------------------------------------------------------------
// CalcExpectedChiSquare()
//--------------------------------------------------------------------------
/**
* <p>Calculates the expected chisq, if applicable.
* <p>Calculates the expected chisq, expected chisq per run, and chisq per run, if applicable.
*
* \param par
* \param totalExpectedChisq expected chisq for all run blocks
@ -130,6 +111,8 @@ void PFitterFcn::CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_
// only do something for chisq
if (fUseChi2) {
Double_t value = 0.0;
// single histo
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i'
expectedChisqPerRun.push_back(value);

98
src/classes/PMsrHandler.cpp
View File

@ -1019,15 +1019,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf",
i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]);
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 (fWriteExpectedChisq)
fout << str.Data() << endl;
@ -1035,6 +1028,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
cout << str.Data() << endl;
}
}
} else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written
for (UInt_t i=0; i<fStatistic.fNdfPerHisto.size(); i++) {
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
}
}
} else { // maxLH
str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
@ -1067,15 +1070,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf",
i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]);
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 (fWriteExpectedChisq)
fout << str.Data() << endl;
@ -1083,6 +1079,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
cout << str.Data() << endl;
}
}
} else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written
for (UInt_t i=0; i<fStatistic.fNdfPerHisto.size(); i++) {
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (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);
@ -1137,15 +1143,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf",
i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]);
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 (fWriteExpectedChisq)
fout << str.Data() << endl;
@ -1153,6 +1152,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
cout << str.Data() << endl;
}
}
} else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written
for (UInt_t i=0; i<fStatistic.fNdfPerHisto.size(); i++) {
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
}
}
} else { // maxLH
str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
@ -1191,15 +1200,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf",
i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]);
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 (fWriteExpectedChisq)
fout << str.Data() << endl;
@ -1207,6 +1209,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
cout << str.Data() << endl;
}
}
} else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written
for (UInt_t i=0; i<fStatistic.fNdfPerHisto.size(); i++) {
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (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);
@ -3970,7 +3982,7 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines)
// GetNoOfFitParameters (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate the number of fit parameters for single histo fits.
* <p>Calculate the number of fit parameters.
*
* \param idx run block index
*/
@ -3993,7 +4005,7 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx)
return 0;
}
// get N0 parameter, possible parameter number or function
// get N0 parameter, possible parameter number or function (single histo fit)
if (fRuns[idx].GetNormParamNo() != -1) {
if (fRuns[idx].GetNormParamNo() < MSR_PARAM_FUN_OFFSET) // parameter
paramVector.push_back(fRuns[idx].GetNormParamNo());
@ -4001,10 +4013,18 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx)
funVector.push_back(fRuns[idx].GetNormParamNo() - MSR_PARAM_FUN_OFFSET);
}
// get background parameter, for the case the background is fitted.
// get background parameter, for the case the background is fitted (single histo fit)
if (fRuns[idx].GetBkgFitParamNo() != -1)
paramVector.push_back(fRuns[idx].GetBkgFitParamNo());
// get alpha parameter if present (asymmetry fit)
if (fRuns[idx].GetAlphaParamNo() != -1)
paramVector.push_back(fRuns[idx].GetAlphaParamNo());
// get beta parameter if present (asymmetry fit)
if (fRuns[idx].GetBetaParamNo() != -1)
paramVector.push_back(fRuns[idx].GetBetaParamNo());
// go through the theory block and collect parameters
// possible entries: number -> parameter, fun<number> -> function, map<number> -> maps
for (UInt_t i=0; i<fTheory.size(); i++) {

23
src/classes/PRunAsymmetry.cpp
View File

@ -219,6 +219,22 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
return chisq;
}
//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate expected chi-square. Currently not implemented since not clear what to be done.
*
* <b>return:</b>
* - chisq value == 0.0
*
* \param par parameter vector iterated by minuit2
*/
Double_t PRunAsymmetry::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
return 0.0;
}
//--------------------------------------------------------------------------
// CalcMaxLikelihood
//--------------------------------------------------------------------------
@ -244,13 +260,6 @@ Double_t PRunAsymmetry::CalcMaxLikelihood(const std::vector<Double_t>& par)
*/
UInt_t PRunAsymmetry::GetNoOfFitBins()
{
// Double_t time;
// fNoOfFitBins=0;
// for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
// time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep();
// if ((time >= fFitStartTime) && (time <= fFitEndTime))
// fNoOfFitBins++;
// }
CalcNoOfFitBins();
return fNoOfFitBins;

199
src/classes/PRunListCollection.cpp
View File

@ -80,7 +80,7 @@ PRunListCollection::~PRunListCollection()
}
//--------------------------------------------------------------------------
// Add
// Add (public)
//--------------------------------------------------------------------------
/**
* <p>Adds a processed set of data to the handler.
@ -150,7 +150,7 @@ void PRunListCollection::SetFitRange(const PDoublePairVector fitRange)
}
//--------------------------------------------------------------------------
// GetSingleHistoChisq
// GetSingleHistoChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of <em>all</em> single histogram runs of a msr-file.
@ -171,29 +171,7 @@ Double_t PRunListCollection::GetSingleHistoChisq(const std::vector<Double_t>& pa
}
//--------------------------------------------------------------------------
// GetSingleHistoChisqExpected
//--------------------------------------------------------------------------
/**
* <p>Calculates expected chi-square of the single histogram with run block index idx of a msr-file.
*
* <b>return:</b>
* - expected chi-square of for a single histogram
*
* \param par fit parameter vector
* \param idx run block index
*/
Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector<Double_t>& par, const UInt_t idx) const
{
Double_t expectedChisq = 0.0;
if (idx < fRunSingleHistoList.size())
expectedChisq = fRunSingleHistoList[idx]->CalcChiSquareExpected(par);
return expectedChisq;
}
//--------------------------------------------------------------------------
// GetAsymmetryChisq
// GetAsymmetryChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of <em>all</em> asymmetry runs of a msr-file.
@ -214,7 +192,7 @@ Double_t PRunListCollection::GetAsymmetryChisq(const std::vector<Double_t>& par)
}
//--------------------------------------------------------------------------
// GetMuMinusChisq
// GetMuMinusChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of <em>all</em> mu minus runs of a msr-file.
@ -235,7 +213,7 @@ Double_t PRunListCollection::GetMuMinusChisq(const std::vector<Double_t>& par) c
}
//--------------------------------------------------------------------------
// GetNonMusrChisq
// GetNonMusrChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of <em>all</em> non-muSR runs of a msr-file.
@ -256,7 +234,109 @@ Double_t PRunListCollection::GetNonMusrChisq(const std::vector<Double_t>& par) c
}
//--------------------------------------------------------------------------
// GetSingleHistoMaximumLikelihood
// GetSingleHistoChisqExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates expected chi-square of the single histogram with run block index idx of a msr-file.
*
* <b>return:</b>
* - expected chi-square of for a single histogram
*
* \param par fit parameter vector
* \param idx run block index
*/
Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector<Double_t>& par, const UInt_t idx) const
{
Double_t expectedChisq = 0.0;
if (idx > fMsrInfo->GetMsrRunList()->size()) {
cerr << ">> PRunListCollection::GetSingleHistoChisqExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl;
return expectedChisq;
}
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
// count how many entries of this fit-type are present up to idx
UInt_t subIdx = 0;
for (UInt_t i=0; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
// return the chisq of the single run
switch (type) {
case PRUN_SINGLE_HISTO:
expectedChisq = fRunSingleHistoList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_ASYMMETRY:
expectedChisq = fRunAsymmetryList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_MU_MINUS:
expectedChisq = fRunMuMinusList[subIdx]->CalcChiSquareExpected(par);
break;
case PRUN_NON_MUSR:
expectedChisq = fRunNonMusrList[subIdx]->CalcChiSquareExpected(par);
break;
default:
break;
}
return expectedChisq;
}
//--------------------------------------------------------------------------
// GetSingleRunChisq (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates chi-square of a single run-block entry of the msr-file.
*
* <b>return:</b>
* - chi-square of single run-block entry with index idx
*
* \param par fit parameter vector
* \param idx run block index
*/
Double_t PRunListCollection::GetSingleRunChisq(const std::vector<Double_t>& par, const UInt_t idx) const
{
Double_t chisq = 0.0;
if (idx > fMsrInfo->GetMsrRunList()->size()) {
cerr << ">> PRunListCollection::GetSingleRunChisq() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl;
return chisq;
}
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
// count how many entries of this fit-type are present up to idx
UInt_t subIdx = 0;
for (UInt_t i=0; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
// return the chisq of the single run
switch (type) {
case PRUN_SINGLE_HISTO:
chisq = fRunSingleHistoList[subIdx]->CalcChiSquare(par);
break;
case PRUN_ASYMMETRY:
chisq = fRunAsymmetryList[subIdx]->CalcChiSquare(par);
break;
case PRUN_MU_MINUS:
chisq = fRunMuMinusList[subIdx]->CalcChiSquare(par);
break;
case PRUN_NON_MUSR:
chisq = fRunNonMusrList[subIdx]->CalcChiSquare(par);
break;
default:
break;
}
return chisq;
}
//--------------------------------------------------------------------------
// GetSingleHistoMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates log max-likelihood of <em>all</em> single histogram runs of a msr-file.
@ -277,7 +357,7 @@ Double_t PRunListCollection::GetSingleHistoMaximumLikelihood(const std::vector<D
}
//--------------------------------------------------------------------------
// GetAsymmetryMaximumLikelihood
// GetAsymmetryMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p> Since it is not clear yet how to handle asymmetry fits with max likelihood
@ -299,7 +379,7 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector<Dou
}
//--------------------------------------------------------------------------
// GetMuMinusMaximumLikelihood
// GetMuMinusMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates log max-likelihood of <em>all</em> mu minus runs of a msr-file.
@ -320,7 +400,7 @@ Double_t PRunListCollection::GetMuMinusMaximumLikelihood(const std::vector<Doubl
}
//--------------------------------------------------------------------------
// GetNonMusrMaximumLikelihood
// GetNonMusrMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p> Since it is not clear yet how to handle non musr fits with max likelihood
@ -342,7 +422,7 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vector<Doubl
}
//--------------------------------------------------------------------------
// GetNoOfBinsFitted
// GetNoOfBinsFitted (public)
//--------------------------------------------------------------------------
/**
* <p>Number of bins in run block idx to be fitted. Only used for single histogram
@ -357,15 +437,44 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const
{
UInt_t result = 0;
if (idx < fRunSingleHistoList.size())
result = fRunSingleHistoList[idx]->GetNoOfFitBins();
if (idx > fMsrInfo->GetMsrRunList()->size()) {
cerr << ">> PRunListCollection::GetNoOfBinsFitted() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl;
return result;
}
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
// count how many entries of this fit-type are present up to idx
UInt_t subIdx = 0;
for (UInt_t i=0; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
// return the chisq of the single run
switch (type) {
case PRUN_SINGLE_HISTO:
result = fRunSingleHistoList[subIdx]->GetNoOfFitBins();
break;
case PRUN_ASYMMETRY:
result = fRunAsymmetryList[subIdx]->GetNoOfFitBins();
break;
case PRUN_MU_MINUS:
result = fRunMuMinusList[subIdx]->GetNoOfFitBins();
break;
case PRUN_NON_MUSR:
result = fRunNonMusrList[subIdx]->GetNoOfFitBins();
break;
default:
break;
}
return result;
}
//--------------------------------------------------------------------------
// GetTotalNoOfBinsFitted
// GetTotalNoOfBinsFitted (public)
//--------------------------------------------------------------------------
/**
* <p>Counts the total number of bins to be fitted.
@ -393,7 +502,7 @@ UInt_t PRunListCollection::GetTotalNoOfBinsFitted() const
}
//--------------------------------------------------------------------------
// GetSingleHisto
// GetSingleHisto (public)
//--------------------------------------------------------------------------
/**
* <p>Get a processed single histogram data set.
@ -436,7 +545,7 @@ PRunData* PRunListCollection::GetSingleHisto(UInt_t index, EDataSwitch tag)
}
//--------------------------------------------------------------------------
// GetAsymmetry
// GetAsymmetry (public)
//--------------------------------------------------------------------------
/**
* <p>Get a processed asymmetry data set.
@ -479,7 +588,7 @@ PRunData* PRunListCollection::GetAsymmetry(UInt_t index, EDataSwitch tag)
}
//--------------------------------------------------------------------------
// GetMuMinus
// GetMuMinus (public)
//--------------------------------------------------------------------------
/**
* <p>Get a processed mu minus data set.
@ -519,7 +628,7 @@ PRunData* PRunListCollection::GetMuMinus(UInt_t index, EDataSwitch tag)
}
//--------------------------------------------------------------------------
// GetNonMusr
// GetNonMusr (public)
//--------------------------------------------------------------------------
/**
* <p>Get a processed non-muSR data set.
@ -559,7 +668,7 @@ PRunData* PRunListCollection::GetNonMusr(UInt_t index, EDataSwitch tag)
}
//--------------------------------------------------------------------------
// GetTemp
// GetTemp (public)
//--------------------------------------------------------------------------
/**
* <p>Get the temperature from the data set.
@ -575,7 +684,7 @@ const PDoublePairVector* PRunListCollection::GetTemp(const TString &runName) con
}
//--------------------------------------------------------------------------
// GetField
// GetField (public)
//--------------------------------------------------------------------------
/**
* <p>Get the magnetic field from the data set.
@ -591,7 +700,7 @@ Double_t PRunListCollection::GetField(const TString &runName) const
}
//--------------------------------------------------------------------------
// GetEnergy
// GetEnergy (public)
//--------------------------------------------------------------------------
/**
* <p>Get the muon implantation energy from the data set.
@ -607,7 +716,7 @@ Double_t PRunListCollection::GetEnergy(const TString &runName) const
}
//--------------------------------------------------------------------------
// GetSetup
// GetSetup (public)
//--------------------------------------------------------------------------
/**
* <p>Get the setup information from the data set.
@ -623,7 +732,7 @@ const Char_t* PRunListCollection::GetSetup(const TString &runName) const
}
//--------------------------------------------------------------------------
// GetXAxisTitle
// GetXAxisTitle (public)
//--------------------------------------------------------------------------
/**
* <p>Get the x-axis title (used with non-muSR fit).
@ -656,7 +765,7 @@ const Char_t* PRunListCollection::GetXAxisTitle(const TString &runName, const UI
}
//--------------------------------------------------------------------------
// GetYAxisTitle
// GetYAxisTitle (public)
//--------------------------------------------------------------------------
/**
* <p>Get the y-axis title (used with non-muSR fit).

16
src/classes/PRunMuMinus.cpp
View File

@ -96,6 +96,22 @@ Double_t PRunMuMinus::CalcChiSquare(const std::vector<Double_t>& par)
return chisq;
}
//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate expected chi-square. Currently not implemented.
*
* <b>return:</b>
* - chisq value == 0.0
*
* \param par parameter vector iterated by minuit2
*/
Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
return 0.0;
}
//--------------------------------------------------------------------------
// CalcMaxLikelihood
//--------------------------------------------------------------------------

16
src/classes/PRunNonMusr.cpp
View File

@ -118,6 +118,22 @@ Double_t PRunNonMusr::CalcChiSquare(const std::vector<Double_t>& par)
return chisq;
}
//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate expected chi-square. Currently not implemented since not clear what to be done.
*
* <b>return:</b>
* - chisq value == 0.0
*
* \param par parameter vector iterated by minuit2
*/
Double_t PRunNonMusr::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
return 0.0;
}
//--------------------------------------------------------------------------
// CalcMaxLikelihood
//--------------------------------------------------------------------------

10
src/classes/PStartupHandler.cpp
View File

@ -166,7 +166,7 @@ void PStartupHandler::OnStartDocument()
fFourierDefaults.fPlotRange[1] = -1.0;
fFourierDefaults.fPhaseIncrement = 1.0;
fWriteExpectedChisq = false;
fWritePerRunBlockChisq = false;
}
//--------------------------------------------------------------------------
@ -195,8 +195,8 @@ void PStartupHandler::OnStartElement(const Char_t *str, const TList *attributes)
{
if (!strcmp(str, "data_path")) {
fKey = eDataPath;
} else if (!strcmp(str, "write_expected_chisq")) {
fKey = eWriteExpectedChisq;
} else if (!strcmp(str, "write_per_run_block_chisq")) {
fKey = eWritePerRunBlockChisq;
} else if (!strcmp(str, "marker")) {
fKey = eMarker;
} else if (!strcmp(str, "color")) {
@ -251,10 +251,10 @@ void PStartupHandler::OnCharacters(const Char_t *str)
// add str to the path list
fDataPathList.push_back(str);
break;
case eWriteExpectedChisq:
case eWritePerRunBlockChisq:
tstr = TString(str);
if (tstr.BeginsWith("y") || tstr.BeginsWith("Y"))
fWriteExpectedChisq = true;
fWritePerRunBlockChisq = true;
break;
case eMarker:
// check that str is a number

2
src/include/PFitterFcn.h
View File

@ -51,7 +51,7 @@ class PFitterFcn : public ROOT::Minuit2::FCNBase
Double_t operator()(const std::vector<Double_t> &par) const;
UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); }
UInt_t GetNoOfFittedBins(const UInt_t idx);
UInt_t GetNoOfFittedBins(const UInt_t idx) { return fRunListCollection->GetNoOfBinsFitted(idx); }
void CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &expectedChisqPerRun);
private:

3
src/include/PMusr.h
View File

@ -583,7 +583,8 @@ typedef struct {
PMsrLines fStatLines; ///< statistics block in msr-file clear text
TString fDate; ///< string holding fitting date and time
Bool_t fChisq; ///< flag telling if min = chi2 or min = max.likelihood
Double_t fMin; ///< chi2 or max. likelihood
Double_t fMin; ///< chisq or max. likelihood
PDoubleVector fMinPerHisto; ///< chisq or max. likelihood per histo
UInt_t fNdf; ///< number of degrees of freedom
Double_t fMinExpected; ///< expected total chi2 or max. likelihood
PDoubleVector fMinExpectedPerHisto; ///< expected pre histo chi2 or max. likelihood

1
src/include/PRunAsymmetry.h
View File

@ -46,6 +46,7 @@ class PRunAsymmetry : public PRunBase
virtual ~PRunAsymmetry();
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory();

4
src/include/PRunListCollection.h
View File

@ -59,11 +59,13 @@ class PRunListCollection
virtual void SetFitRange(const PDoublePairVector fitRange);
virtual Double_t GetSingleHistoChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetSingleHistoChisqExpected(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetAsymmetryChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetMuMinusChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetNonMusrChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetSingleHistoChisqExpected(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetSingleRunChisq(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetSingleHistoMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetAsymmetryMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetMuMinusMaximumLikelihood(const std::vector<Double_t>& par) const;

1
src/include/PRunMuMinus.h
View File

@ -45,6 +45,7 @@ class PRunMuMinus : public PRunBase
virtual ~PRunMuMinus();
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory();

1
src/include/PRunNonMusr.h
View File

@ -46,6 +46,7 @@ class PRunNonMusr : public PRunBase
virtual ~PRunNonMusr();
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory();

6
src/include/PStartupHandler.h
View File

@ -75,14 +75,14 @@ class PStartupHandler : public TObject, public TQObject
virtual void CheckLists();
virtual const Bool_t GetWriteExpectedChisq() { return fWriteExpectedChisq; } ///< returns the write_expected_chisq flag
virtual const Bool_t GetWritePerRunBlockChisq() { return fWritePerRunBlockChisq; } ///< returns the write_per_run_block_chisq flag
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
private:
enum EKeyWords {eEmpty, eComment, eDataPath, eWriteExpectedChisq,
enum EKeyWords {eEmpty, eComment, eDataPath, eWritePerRunBlockChisq,
eFourierSettings, eUnits, eFourierPower, eApodization, ePlot, ePhase, ePhaseIncrement,
eRootSettings, eMarkerList, eMarker,
eColorList, eColor};
@ -90,7 +90,7 @@ class PStartupHandler : public TObject, public TQObject
Bool_t fStartupFileFound; ///< startup file found flag
TString fStartupFilePath; ///< full musrfit_startup.xml startup file paths
Bool_t fWriteExpectedChisq; ///< flag showing if the expected chisq shall be written to the msr-file
Bool_t fWritePerRunBlockChisq; ///< flag showing if per run block chisq and the expected chisq shall be written to the msr-file
PMsrFourierStructure fFourierDefaults; ///< Fourier defaults
PStringVector fDataPathList; ///< search data path list
PIntVector fMarkerList; ///< marker list

2
src/musrfit.cpp
View File

@ -463,7 +463,7 @@ int main(int argc, char *argv[])
}
// read msr-file
PMsrHandler *msrHandler = new PMsrHandler(filename, startupHandler->GetWriteExpectedChisq());
PMsrHandler *msrHandler = new PMsrHandler(filename, startupHandler->GetWritePerRunBlockChisq());
status = msrHandler->ReadMsrFile();
if (status != PMUSR_SUCCESS) {
switch (status) {

2
src/musrfit_startup.xml
View File

@ -13,7 +13,7 @@
<data_path>/afs/psi.ch/project/bulkmusr/data/gpd</data_path>
<data_path>/afs/psi.ch/project/bulkmusr/data/ltf</data_path>
<data_path>/afs/psi.ch/project/bulkmusr/data/alc</data_path>
<write_expected_chisq>n</write_expected_chisq>
<write_per_run_block_chisq>n</write_per_run_block_chisq>
<fourier_settings>
<units>Gauss</units>
<fourier_power>0</fourier_power>