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

added expected chisq calculation for single histogram fits (MUSR-194)

Browse Source
This commit is contained in:
nemu 2011-07-09 16:26:42 +00:00
parent 0ab456e5f5
commit d9f1f1235e
18 changed files with 606 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ChangeLog
View File

@ -6,6 +6,9 @@
changes since 0.9.0
===================================
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)
NEW the chi^2 calculation in single-histogram and asymmetry fits is parallelized
if musrfit is built using a compiler supporting OpenMP (e.g. GCC >= 4.2)
Using --disable-omp this feature can be disabled on the configure level.

53
src/classes/PFitter.cpp
View File

@ -170,7 +170,27 @@ Bool_t PFitter::DoFit()
Int_t ndf = fFitterFcn->GetTotalNoOfFittedBins() - usedParams;
Double_t val = (*fFitterFcn)(param);
if (fUseChi2) {
// calculate expected chisq
Double_t totalExpectedChisq = 0.0;
std::vector<Double_t> expectedChisqPerHisto;
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
cout << endl << endl << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf;
if (totalExpectedChisq != 0.0) {
cout << endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf;
UInt_t ndf_histo = 0;
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];
}
// clean up
expectedChisqPerHisto.clear();
}
} else { // max. log likelihood
cout << endl << endl << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf;
}
@ -1421,6 +1441,39 @@ Bool_t PFitter::ExecuteSave()
return false;
}
// handle expected chisq if applicable
if (fUseChi2) {
fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back
// calculate expected chisq
std::vector<Double_t> param;
Double_t totalExpectedChisq = 0.0;
std::vector<Double_t> expectedChisqPerHisto;
std::vector<UInt_t> ndfPerHisto;
for (UInt_t i=0; i<fParams.size(); i++)
param.push_back(fParams[i].fValue);
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only
// get the ndf's of the histos
UInt_t ndf_histo;
for (UInt_t i=0; i<expectedChisqPerHisto.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->fMinExpected = totalExpectedChisq;
statistics->fMinExpectedPerHisto = expectedChisqPerHisto;
statistics->fNdfPerHisto = ndfPerHisto;
}
}
}
cout << ">> PFitter::ExecuteSave(): will write minuit2 output file ..." << endl;
ofstream fout;

47
src/classes/PFitterFcn.cpp
View File

@ -65,6 +65,8 @@ PFitterFcn::~PFitterFcn()
{
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Minuit2 interface function call routine. This is the function which should be minimized.
@ -90,3 +92,48 @@ 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.
*
* \param par
* \param totalExpectedChisq expected chisq for all run blocks
* \param expectedChisqPerRun expected chisq vector for all the run blocks
*/
void PFitterFcn::CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &expectedChisqPerRun)
{
// init expected chisq related variables
totalExpectedChisq = 0.0;
expectedChisqPerRun.clear();
// only do something for chisq
if (fUseChi2) {
Double_t value = 0.0;
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);
totalExpectedChisq += value;
}
}
}

392
src/classes/PMsrHandler.cpp
View File

@ -52,7 +52,8 @@ using namespace std;
*
* \param fileName name of a msr-file.
*/
PMsrHandler::PMsrHandler(const Char_t *fileName) : fFileName(fileName)
PMsrHandler::PMsrHandler(const Char_t *fileName, const Bool_t writeExpectedChisq) :
fWriteExpectedChisq(writeExpectedChisq), fFileName(fileName)
{
// init variables
fMsrBlockCounter = 0;
@ -64,6 +65,9 @@ PMsrHandler::PMsrHandler(const Char_t *fileName) : fFileName(fileName)
fStatistic.fChisq = true;
fStatistic.fMin = -1.0;
fStatistic.fNdf = 0;
fStatistic.fMinExpected = 0.0;
fStatistic.fMinExpectedPerHisto.clear();
fStatistic.fNdfPerHisto.clear();
fFuncHandler = 0;
@ -95,6 +99,8 @@ PMsrHandler::~PMsrHandler()
fCommands.clear();
fPlots.clear();
fStatistic.fStatLines.clear();
fStatistic.fMinExpectedPerHisto.clear();
fStatistic.fNdfPerHisto.clear();
fParamInUse.clear();
if (fFuncHandler) {
@ -996,21 +1002,46 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
} else if (sstr.BeginsWith("chisq") || sstr.BeginsWith("maxLH")) {
partialStatisticBlockFound = false;
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq)
str = " chisq = ";
else
str = " maxLH = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
if (fStatistic.fChisq)
str += ", chisq/NDF = ";
else
str += ", maxLH/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
if (fStatistic.fChisq) {
str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
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]);
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);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
}
} else {
fout << "*** FIT DID NOT CONVERGE ***" << endl;
if (messages)
@ -1020,22 +1051,41 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
partialStatisticBlockFound = false;
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq) { // chisq
str = " chisq = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
str += ", chisq/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
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]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
}
}
}
} else { // max. log. liklihood
str = " maxLH = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
str += ", maxLH/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
@ -1046,11 +1096,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
cout << endl << "*** FIT DID NOT CONVERGE ***" << endl;
}
} else {
if (str.Length() > 0)
fout << str.Data() << endl;
else // only write endl if not eof is reached. This is preventing growing msr-files, i.e. more and more empty lines at the end of the file
if (str.Length() > 0) {
sstr = str;
sstr.Remove(TString::kLeading, ' ');
if (!sstr.BeginsWith("expected chisq") && !sstr.BeginsWith("run block"))
fout << str.Data() << endl;
} else { // only write endl if not eof is reached. This is preventing growing msr-files, i.e. more and more empty lines at the end of the file
if (!fin.eof())
fout << endl;
}
}
break;
default:
@ -1066,28 +1120,47 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
TDatime dt;
fout << "STATISTIC --- " << dt.AsSQLString() << endl;
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq) { // chisq
str = " chisq = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
str += ", chisq/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
if (fStatistic.fChisq) {
str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
} else { // max. log. liklihood
str = " maxLH = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
str += ", maxLH/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
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]);
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);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
}
} else {
} else {
fout << "*** FIT DID NOT CONVERGE ***" << endl;
if (messages)
cout << endl << "*** FIT DID NOT CONVERGE ***" << endl;
@ -1102,28 +1175,47 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout << "STATISTIC --- " << dt.AsSQLString() << endl;
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq) { // chisq
str = " chisq = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
str += ", chisq/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
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]);
if (fWriteExpectedChisq)
fout << str.Data() << endl;
if (messages)
cout << str.Data() << endl;
}
}
}
} else { // max. log. liklihood
str = " maxLH = ";
str += fStatistic.fMin;
str += ", NDF = ";
str += fStatistic.fNdf;
str += ", maxLH/NDF = ";
str += fStatistic.fMin / fStatistic.fNdf;
str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
if (messages)
cout << endl << str.Data() << endl;
}
} else {
fout << "*** FIT DID NOT CONVERGE ***" << endl;
fout << "*** FIT DID NOT CONVERGE (4) ***" << endl;
if (messages)
cout << endl << "*** FIT DID NOT CONVERGE ***" << endl;
}
@ -3813,7 +3905,8 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines)
tstr.Remove(TString::kLeading, ' ');
if (tstr.Length() > 0) {
if (!tstr.BeginsWith("#") && !tstr.BeginsWith("STATISTIC") && !tstr.BeginsWith("chisq") &&
!tstr.BeginsWith("maxLH") && !tstr.BeginsWith("*** FIT DID NOT CONVERGE ***")) {
!tstr.BeginsWith("maxLH") && !tstr.BeginsWith("*** FIT DID NOT CONVERGE ***") &&
!tstr.BeginsWith("expected chisq") && !tstr.BeginsWith("run block")) {
cerr << endl << ">> PMsrHandler::HandleStatisticEntry: **SYNTAX ERROR** in line " << lines[i].fLineNo;
cerr << endl << ">> '" << lines[i].fLine.Data() << "'";
cerr << endl << ">> not a valid STATISTIC block line";
@ -3872,6 +3965,185 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines)
return true;
}
//--------------------------------------------------------------------------
// GetNoOfFitParameters (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate the number of fit parameters for single histo fits.
*
* \param idx run block index
*/
UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx)
{
UInt_t noOfFitParameters = 0;
PUIntVector paramVector;
PUIntVector funVector;
PUIntVector mapVector;
TObjArray *tokens = 0;
TObjString *ostr = 0;
TString str;
UInt_t k, dval;
Int_t status, pos;
// check that idx is valid
if (idx >= fRuns.size()) {
cerr << endl << ">> PMsrHandler::GetNoOfFitParameters() **ERROR** idx=" << idx << ", out of range fRuns.size()=" << fRuns.size();
cerr << endl;
return 0;
}
// get N0 parameter, possible parameter number or function
if (fRuns[idx].GetNormParamNo() != -1) {
if (fRuns[idx].GetNormParamNo() < MSR_PARAM_FUN_OFFSET) // parameter
paramVector.push_back(fRuns[idx].GetNormParamNo());
else // function
funVector.push_back(fRuns[idx].GetNormParamNo() - MSR_PARAM_FUN_OFFSET);
}
// get background parameter, for the case the background is fitted.
if (fRuns[idx].GetBkgFitParamNo() != -1)
paramVector.push_back(fRuns[idx].GetBkgFitParamNo());
// 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++) {
// remove potential comments
str = fTheory[i].fLine;
pos = str.Index('#');
if (pos >= 0)
str.Resize(pos);
// tokenize
tokens = str.Tokenize(" \t");
if (!tokens) {
mapVector.clear();
funVector.clear();
paramVector.clear();
return 0;
}
for (Int_t j=0; j<tokens->GetEntries(); j++) {
ostr = dynamic_cast<TObjString*>(tokens->At(j));
str = ostr->GetString();
// check for parameter number
if (str.IsDigit()) {
dval = str.Atoi();
paramVector.push_back(dval);
}
// check for map
if (str.Contains("map")) {
status = sscanf(str.Data(), "map%d", &dval);
if (status == 1) {
mapVector.push_back(dval);
}
}
// check for function
if (str.Contains("fun")) {
status = sscanf(str.Data(), "fun%d", &dval);
if (status == 1) {
funVector.push_back(dval);
}
}
}
delete tokens;
tokens = 0;
}
// go through the function block and collect parameters
for (UInt_t i=0; i<funVector.size(); i++) {
// find the proper function in the function block
for (k=0; k<fFunctions.size(); k++) {
status = sscanf(fFunctions[k].fLine.Data(), "fun%d", &dval);
if (status == 1) {
if (dval == funVector[i])
break;
}
}
// check if everything has been found at all
if (k == fFunctions.size()) {
cerr << endl << ">> PMsrHandler::GetNoOfFitParameters() **ERROR** couldn't find fun" << funVector[i];
cerr << endl << endl;
// clean up
mapVector.clear();
funVector.clear();
paramVector.clear();
return 0;
}
// remove potential comments
str = fFunctions[k].fLine;
pos = str.Index('#');
if (pos >= 0)
str.Resize(pos);
// tokenize
tokens = str.Tokenize(" \t");
if (!tokens) {
mapVector.clear();
funVector.clear();
paramVector.clear();
return 0;
}
// filter out parameters and maps
for (Int_t j=0; j<tokens->GetEntries(); j++) {
ostr = dynamic_cast<TObjString*>(tokens->At(j));
str = ostr->GetString();
// check for parameter
if (str.BeginsWith("par")) {
status = sscanf(str.Data(), "par%d", &dval);
if (status == 1)
paramVector.push_back(dval);
}
// check for map
if (str.BeginsWith("map")) {
status = sscanf(str.Data(), "map%d", &dval);
if (status == 1)
mapVector.push_back(dval);
}
}
}
// go through the map and collect parameters
for (UInt_t i=0; i<mapVector.size(); i++) {
paramVector.push_back(fRuns[idx].GetMap(mapVector[i]-1));
}
// eliminated multiple identical entries in paramVector
PUIntVector param;
param.push_back(paramVector[0]);
for (UInt_t i=0; i<paramVector.size(); i++) {
for (k=0; k<param.size(); k++) {
if (param[k] == paramVector[i])
break;
}
if (k == param.size())
param.push_back(paramVector[i]);
}
// calculate the number of fit parameters with step != 0
for (UInt_t i=0; i<param.size(); i++) {
if (fParam[param[i]-1].fStep != 0.0)
noOfFitParameters++;
}
// cleanup
param.clear();
mapVector.clear();
funVector.clear();
paramVector.clear();
return noOfFitParameters;
}
//--------------------------------------------------------------------------
// FillParameterInUse (private)
//--------------------------------------------------------------------------

2
src/classes/PRunBase.cpp
View File

@ -123,8 +123,6 @@ PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo,
*/
PRunBase::~PRunBase()
{
fParamNo.clear();
fT0s.clear();
fFuncValues.clear();

45
src/classes/PRunListCollection.cpp
View File

@ -170,6 +170,28 @@ Double_t PRunListCollection::GetSingleHistoChisq(const std::vector<Double_t>& pa
return chisq;
}
//--------------------------------------------------------------------------
// 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
//--------------------------------------------------------------------------
@ -319,6 +341,29 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vector<Doubl
return mlh;
}
//--------------------------------------------------------------------------
// GetNoOfBinsFitted
//--------------------------------------------------------------------------
/**
* <p>Number of bins in run block idx to be fitted. Only used for single histogram
* fitting together with the expected chisq.
*
* <b>return:</b>
* - number of bins fitted.
*
* \param idx run block index
*/
UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const
{
UInt_t result = 0;
if (idx < fRunSingleHistoList.size())
result = fRunSingleHistoList[idx]->GetNoOfFitBins();
return result;
}
//--------------------------------------------------------------------------
// GetTotalNoOfBinsFitted
//--------------------------------------------------------------------------

101
src/classes/PRunSingleHisto.cpp
View File

@ -185,6 +185,99 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
return chisq;
}
//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate expected chi-square.
*
* <b>return:</b>
* - chisq value
*
* \param par parameter vector iterated by minuit2
*/
Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
Double_t chisq = 0.0;
Double_t diff = 0.0;
Double_t theo = 0.0;
Double_t N0 = 0.0;
// check if norm is a parameter or a function
if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
N0 = par[fRunInfo->GetNormParamNo()-1];
} else { // norm is a function
// get function number
UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par);
}
// get tau
Double_t tau;
if (fRunInfo->GetLifetimeParamNo() != -1)
tau = par[fRunInfo->GetLifetimeParamNo()-1];
else
tau = PMUON_LIFETIME;
// get background
Double_t bkg;
if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
bkg = fBackground;
} else { // fixed bkg given
bkg = fRunInfo->GetBkgFix(0);
}
} else { // bkg fitted
bkg = par[fRunInfo->GetBkgFitParamNo()-1];
}
// calculate functions
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
Int_t funcNo = fMsrInfo->GetFuncNo(i);
fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par);
}
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
// for a given set of parameters---which should be done outside of the parallelized loop.
// For all other functions it means a tiny and acceptable overhead.
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;
diff = fData.GetValue()->at(i) - theo;
chisq += diff*diff / theo;
}
// the correction factor is need since the data scales like pack*t_res,
// whereas the error scales like sqrt(pack*t_res)
if (fScaleN0AndBkg)
chisq *= fRunInfo->GetPacking() * (fTimeResolution * 1.0e3);
return chisq;
}
//--------------------------------------------------------------------------
// CalcMaxLikelihood (public)
//--------------------------------------------------------------------------
@ -360,14 +453,6 @@ void PRunSingleHisto::CalcTheory()
*/
UInt_t PRunSingleHisto::GetNoOfFitBins()
{
// fNoOfFitBins=0;
//
// Double_t time;
// 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;

9
src/classes/PStartupHandler.cpp
View File

@ -165,6 +165,8 @@ void PStartupHandler::OnStartDocument()
fFourierDefaults.fPlotRange[0] = -1.0;
fFourierDefaults.fPlotRange[1] = -1.0;
fFourierDefaults.fPhaseIncrement = 1.0;
fWriteExpectedChisq = false;
}
//--------------------------------------------------------------------------
@ -193,6 +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, "marker")) {
fKey = eMarker;
} else if (!strcmp(str, "color")) {
@ -247,6 +251,11 @@ void PStartupHandler::OnCharacters(const Char_t *str)
// add str to the path list
fDataPathList.push_back(str);
break;
case eWriteExpectedChisq:
tstr = TString(str);
if (tstr.BeginsWith("y") || tstr.BeginsWith("Y"))
fWriteExpectedChisq = true;
break;
case eMarker:
// check that str is a number
tstr = TString(str);

2
src/classes/PTheory.cpp
View File

@ -307,8 +307,6 @@ PTheory::PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent) : f
fUserParam.resize(fParamNo.size());
}
//cout << endl << "debug> fUserFcn=" << fUserFcn << ", fUserFcn->NeedGlobalPart()=" << fUserFcn->NeedGlobalPart() << ", gGlobalUserFcn=" << gGlobalUserFcn << endl;
// check if the global part of the user function is needed
if (fUserFcn->NeedGlobalPart()) {
fUserFcn->SetGlobalPart(gGlobalUserFcn, fUserFcnIdx); // invoke or retrieve global user function object

4
src/include/PFitterFcn.h
View File

@ -48,9 +48,11 @@ class PFitterFcn : public ROOT::Minuit2::FCNBase
~PFitterFcn();
Double_t Up() const { return fUp; }
Double_t operator()(const std::vector<Double_t>&) const;
Double_t operator()(const std::vector<Double_t> &par) const;
UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); }
UInt_t GetNoOfFittedBins(const UInt_t idx);
void CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &expectedChisqPerRun);
private:
Double_t fUp; ///< for chisq == 1.0, i.e. errors are 1 std. deviation errors. for log max-likelihood == 0.5, i.e. errors are 1 std. deviation errors (for details see the minuit2 user manual).

5
src/include/PMsrHandler.h
View File

@ -47,7 +47,7 @@
class PMsrHandler
{
public:
PMsrHandler(const Char_t *fileName);
PMsrHandler(const Char_t *fileName, const Bool_t writeExpectedChisq=false);
virtual ~PMsrHandler();
virtual Int_t ReadMsrFile();
@ -95,6 +95,7 @@ class PMsrHandler
virtual Double_t EvalFunc(UInt_t i, vector<Int_t> map, vector<Double_t> param)
{ return fFuncHandler->Eval(i,map,param); }
virtual UInt_t GetNoOfFitParameters(UInt_t idx);
virtual Int_t ParameterInUse(UInt_t paramNo);
virtual Bool_t CheckRunBlockIntegrity();
virtual Bool_t CheckUniquenessOfParamNames(UInt_t &parX, UInt_t &parY);
@ -105,6 +106,8 @@ class PMsrHandler
virtual void CheckMaxLikelihood();
private:
Bool_t fWriteExpectedChisq; ///< flag shows if expected chisq shall be written to the msr-file
TString fFileName; ///< file name of the msr-file
TString fMsrFileDirectoryPath; ///< msr-file directory path
TString fTitle; ///< holds the title string of the msr-file

7
src/include/PMusr.h
View File

@ -582,9 +582,12 @@ typedef struct {
Bool_t fValid; ///< flag showing if the statistics block is valid, i.e. a fit took place which converged
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.likelyhood
Double_t fMin; ///< chi2 or max. likelyhood
Bool_t fChisq; ///< flag telling if min = chi2 or min = max.likelihood
Double_t fMin; ///< chi2 or max. likelihood
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
PUIntVector fNdfPerHisto; ///< number of degrees of freedom per histo
} PMsrStatisticStructure;
//-------------------------------------------------------------

2
src/include/PRunBase.h
View File

@ -74,8 +74,6 @@ class PRunBase
PMsrRunBlock *fRunInfo; ///< run info used to filter out needed infos of a run
PRunDataHandler *fRawData; ///< holds the raw run data
PIntVector fParamNo; ///< vector of parameter numbers for the specifc run
PRunData fData; ///< data to be fitted, viewed, i.e. binned data
Double_t fTimeResolution; ///< time resolution in (us)
PIntVector fT0s; ///< all t0's of a run! The derived classes will handle it

10
src/include/PRunListCollection.h
View File

@ -59,6 +59,7 @@ 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;
@ -68,6 +69,7 @@ class PRunListCollection
virtual Double_t GetMuMinusMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetNonMusrMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const;
virtual UInt_t GetTotalNoOfBinsFitted() const;
virtual UInt_t GetNoOfSingleHisto() const { return fRunSingleHistoList.size(); } ///< returns the number of single histogram data sets present in the msr-file
@ -91,10 +93,10 @@ class PRunListCollection
PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler
PRunDataHandler *fData; ///< pointer to the run-data handler
vector<PRunSingleHisto*> fRunSingleHistoList; ///< stores all precessed single histogram data
vector<PRunAsymmetry*> fRunAsymmetryList; ///< stores all precessed asymmetry data
vector<PRunMuMinus*> fRunMuMinusList; ///< stores all precessed mu-minus data
vector<PRunNonMusr*> fRunNonMusrList; ///< stores all precessed non-muSR data
vector<PRunSingleHisto*> fRunSingleHistoList; ///< stores all processed single histogram data
vector<PRunAsymmetry*> fRunAsymmetryList; ///< stores all processed asymmetry data
vector<PRunMuMinus*> fRunMuMinusList; ///< stores all processed mu-minus data
vector<PRunNonMusr*> fRunNonMusrList; ///< stores all processed non-muSR data
};
#endif // _PRUNLISTCOLLECTION_H_

1
src/include/PRunSingleHisto.h
View File

@ -45,6 +45,7 @@ class PRunSingleHisto : public PRunBase
virtual ~PRunSingleHisto();
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/PStartupHandler.h
View File

@ -75,13 +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 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,
enum EKeyWords {eEmpty, eComment, eDataPath, eWriteExpectedChisq,
eFourierSettings, eUnits, eFourierPower, eApodization, ePlot, ePhase, ePhaseIncrement,
eRootSettings, eMarkerList, eMarker,
eColorList, eColor};
@ -89,6 +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
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);
PMsrHandler *msrHandler = new PMsrHandler(filename, startupHandler->GetWriteExpectedChisq());
status = msrHandler->ReadMsrFile();
if (status != PMUSR_SUCCESS) {
switch (status) {

1
src/musrfit_startup.xml
View File

@ -13,6 +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>
<fourier_settings>
<units>Gauss</units>
<fourier_power>0</fourier_power>