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

for maxLH it is now possible to get the per-run-block maxLH, and the expected maxLH

Browse Source
This commit is contained in:
suter_a 2016-12-18 10:36:02 +01:00
parent 65d40cfe97
commit c768c27898
9 changed files with 433 additions and 192 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

104
src/classes/PFitter.cpp
View File

@ -1564,64 +1564,66 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave)
}
// handle expected chisq if applicable
if (fUseChi2) {
fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back
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;
// 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);
for (UInt_t i=0; i<fParams.size(); i++)
param.push_back(fParams[i].fValue);
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
// CalcExpectedChiSquare handles both, chisq and mlh
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
// calculate chisq per run
std::vector<Double_t> chisqPerHisto;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
// calculate chisq per run
std::vector<Double_t> chisqPerHisto;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
if (fUseChi2)
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;
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->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();
else
chisqPerHisto.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(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;
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->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;
ofstream fout;

11
src/classes/PFitterFcn.cpp
View File

@ -110,15 +110,20 @@ void PFitterFcn::CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_
totalExpectedChisq = 0.0;
expectedChisqPerRun.clear();
// only do something for chisq
Double_t value = 0.0;
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);
totalExpectedChisq += value;
}
} else { // log max. likelihood
// single histo
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
value = fRunListCollection->GetSingleHistoMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i'
expectedChisqPerRun.push_back(value);
totalExpectedChisq += value;
}
}
}

305
src/classes/PMsrHandler.cpp
View File

@ -825,13 +825,7 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout << left << "lifetime";
fout << fRuns[runNo].GetLifetimeParamNo() << endl;
} else if (sstr.BeginsWith("lifetimecorrection")) {
/* obsolate, hence do nothing here
if ((fRuns[runNo].IsLifetimeCorrected()) &&
((fRuns[runNo].GetFitType() == MSR_FITTYPE_SINGLE_HISTO) ||
(fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO))) {
fout << "lifetimecorrection" << endl;
}
*/
// obsolate, hence do nothing here
} else if (sstr.BeginsWith("map")) {
fout << "map ";
for (UInt_t j=0; j<fRuns[runNo].GetMap()->size(); j++) {
@ -1215,38 +1209,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq) {
str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf);
fout << str.Data() << endl;
} else {
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;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
if (fStatistic.fChisq) {
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
} else {
str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
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) {
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
if (fStatistic.fChisq) {
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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) = (%d/%lf/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
}
} 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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
@ -1256,11 +1250,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t 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 if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written
for (UInt_t i=0; i<fStatistic.fNdfPerHisto.size(); i++) {
if (fStatistic.fChisq) {
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
} else {
str.Form(" run block %d: (NDF/maxLH.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
}
}
} else {
fout << "*** FIT DID NOT CONVERGE ***" << endl;
@ -1272,38 +1278,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq) { // chisq
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 {
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;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
if (fStatistic.fChisq) { // chisq
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
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) {
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
if (fStatistic.fChisq) { // chisq
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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) = (%d/%lf/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
}
} 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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
@ -1313,11 +1319,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t 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);
fout << str.Data() << endl;
if (messages)
cout << endl << 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++) {
if (fStatistic.fChisq) { // chisq
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
} else {
str.Form(" run block %d: (NDF/red.maxLH) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
}
}
} else {
fout << "*** FIT DID NOT CONVERGE ***" << endl;
@ -1328,7 +1346,7 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
if (str.Length() > 0) {
sstr = str;
sstr.Remove(TString::kLeading, ' ');
if (!sstr.BeginsWith("expected chisq") && !sstr.BeginsWith("run block"))
if (!sstr.BeginsWith("expected chisq") && !sstr.BeginsWith("expected maxLH") && !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())
@ -1351,38 +1369,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
if (fStatistic.fValid) { // valid fit result
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 {
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;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
if (fStatistic.fChisq) {
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
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) {
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
if (fStatistic.fChisq) {
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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) = (%d/%lf/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
}
} 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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
@ -1392,13 +1410,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t 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 if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written
for (UInt_t i=0; i<fStatistic.fNdfPerHisto.size(); i++) {
if (fStatistic.fChisq) {
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
} else {
str.Form(" run block %d: (NDF/red.maxLH) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
// check if per run block maxLH needs to be written
if (messages)
cout << str.Data() << endl;
}
}
} else {
fout << "*** FIT DID NOT CONVERGE ***" << endl;
@ -1416,38 +1444,38 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
if (fStatistic.fValid) { // valid fit result
if (fStatistic.fChisq) { // chisq
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 {
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;
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
// check if expected chisq needs to be written
if (fStatistic.fMinExpected != 0.0) {
if (fStatistic.fChisq) { // chisq
str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" expected maxLH = %.1lf, NDF = %d, expected maxLH/NDF = %lf",
fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
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) {
for (UInt_t i=0; i<fStatistic.fMinExpectedPerHisto.size(); i++) {
if (fStatistic.fNdfPerHisto[i] > 0) {
if (fStatistic.fChisq) { // chisq
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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
} else {
str.Form(" run block %d: (NDF/red.maxLH/red.maxLH_e) =(%d/%lf/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
}
} 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 (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
@ -1457,11 +1485,23 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t 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);
fout << str.Data() << endl;
if (messages)
cout << endl << 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++) {
if (fStatistic.fChisq) { // chisq
str.Form(" run block %d: (NDF/red.chisq) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
} else {
str.Form(" run block %d: (NDF/red.maxLH) = (%d/%lf)",
i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i]);
}
if (fStartupOptions) {
if (fStartupOptions->writeExpectedChisq)
fout << str.Data() << endl;
}
if (messages)
cout << str.Data() << endl;
}
}
} else {
fout << "*** FIT DID NOT CONVERGE (4) ***" << endl;
@ -4675,7 +4715,8 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines)
if (tstr.Length() > 0) {
if (!tstr.BeginsWith("#") && !tstr.BeginsWith("STATISTIC") && !tstr.BeginsWith("chisq") &&
!tstr.BeginsWith("maxLH") && !tstr.BeginsWith("*** FIT DID NOT CONVERGE ***") &&
!tstr.BeginsWith("expected chisq") && !tstr.BeginsWith("run block")) {
!tstr.BeginsWith("expected chisq") && !tstr.BeginsWith("expected maxLH") &&
!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";
@ -6059,6 +6100,7 @@ Bool_t PMsrHandler::EstimateN0()
/**
* <p>returns alpha to estimate N0
*/
/*as
Double_t PMsrHandler::GetAlphaEstimateN0()
{
if (fStartupOptions == 0)
@ -6066,6 +6108,7 @@ Double_t PMsrHandler::GetAlphaEstimateN0()
return fStartupOptions->alphaEstimateN0;
}
*/
//--------------------------------------------------------------------------
// NeededPrecision (private)

93
src/classes/PRunListCollection.cpp
View File

@ -355,7 +355,7 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector<Doubl
}
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
if (type == -1) { // i.e. not forun in the RUN block, try the GLOBAL block
if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block
type = fMsrInfo->GetMsrGlobal()->GetFitType();
}
@ -583,6 +583,97 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vector<Doubl
return mlh;
}
//--------------------------------------------------------------------------
// GetSingleHistoMaximumLikelihoodExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates expected mlh of the single histogram with run block index idx of a msr-file.
*
* <b>return:</b>
* - expected mlh of for a single histogram
*
* \param par fit parameter vector
* \param idx run block index
*/
Double_t PRunListCollection::GetSingleHistoMaximumLikelihoodExpected(const std::vector<Double_t>& par, const UInt_t idx) const
{
Double_t expected_mlh = 0.0;
if (idx > fMsrInfo->GetMsrRunList()->size()) {
cerr << ">> PRunListCollection::GetSingleHistoMaximumLikelihoodExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl;
return expected_mlh;
}
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block
type = fMsrInfo->GetMsrGlobal()->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 mlh of the single run
switch (type) {
case PRUN_SINGLE_HISTO:
expected_mlh = fRunSingleHistoList[subIdx]->CalcMaxLikelihoodExpected(par);
break;
default:
break;
}
return expected_mlh;
}
//--------------------------------------------------------------------------
// GetSingleRunMaximumLikelihood (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates mlh of a single run-block entry of the msr-file.
*
* <b>return:</b>
* - mlh of single run-block entry with index idx
*
* \param par fit parameter vector
* \param idx run block index
*/
Double_t PRunListCollection::GetSingleRunMaximumLikelihood(const std::vector<Double_t>& par, const UInt_t idx) const
{
Double_t mlh = 0.0;
if (idx > fMsrInfo->GetMsrRunList()->size()) {
cerr << ">> PRunListCollection::GetSingleRunMaximumLikelihood() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl;
return mlh;
}
Int_t subIdx = 0;
Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType();
if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block
type = fMsrInfo->GetMsrGlobal()->GetFitType();
subIdx = idx;
} else { // found in the RUN block
// count how many entries of this fit-type are present up to idx
for (UInt_t i=0; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
}
// return the mlh of the single run
switch (type) {
case PRUN_SINGLE_HISTO:
mlh = fRunSingleHistoList[subIdx]->CalcMaxLikelihood(par);
break;
default:
break;
}
return mlh;
}
//--------------------------------------------------------------------------
// GetNoOfBinsFitted (public)
//--------------------------------------------------------------------------

100
src/classes/PRunSingleHisto.cpp
View File

@ -393,6 +393,103 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
return 2.0*mllh;
}
//--------------------------------------------------------------------------
// CalcMaxLikelihoodExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate expected log maximum-likelihood.
*
* <b>return:</b>
* - log maximum-likelihood value
*
* \param par parameter vector iterated by minuit2
*/
Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector<Double_t>& par)
{
Double_t mllh = 0.0; // maximum log likelihood assuming poisson distribution for the single bin
Double_t N0;
// 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 maximum log likelihood
Double_t theo;
Double_t data;
Double_t time(1.0);
Int_t i;
// norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns
Double_t normalizer = 1.0;
if (fScaleN0AndBkg)
normalizer = fPacking * (fTimeResolution * 1.0e3);
// 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 = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,theo,data) schedule(dynamic,chunk) reduction(-:mllh)
#endif
for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
// calculate theory for the given parameter set
theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;
theo *= normalizer;
data = normalizer*fData.GetValue()->at(i);
if (theo <= 0.0) {
cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << endl;
continue;
}
if (data > 1.0e-9) { // is this correct?? needs to be checked. See G-test
mllh += data*log(data/theo);
}
}
return 2.0*mllh;
}
//--------------------------------------------------------------------------
// CalcTheory (public)
//--------------------------------------------------------------------------
@ -1571,7 +1668,6 @@ void PRunSingleHisto::EstimateN0()
Double_t tau = PMUON_LIFETIME;
UInt_t t0 = (UInt_t)round(fT0s[0]);
Double_t alpha = fMsrInfo->GetAlphaEstimateN0();
Double_t dval = 0.0;
Double_t nom = 0.0;
Double_t denom = 0.0;
@ -1580,14 +1676,12 @@ void PRunSingleHisto::EstimateN0()
// calc nominator
for (UInt_t i=t0; i<fForward.size(); i++) {
xx = exp(-dt*(Double_t)(i-t0)/tau);
xx += alpha;
nom += xx;
}
// calc denominator
for (UInt_t i=t0; i<fForward.size(); i++) {
xx = exp(-dt*(Double_t)(i-t0)/tau);
xx += alpha;
dval = fForward[i];
if (dval > 0)
denom += xx*xx/dval;

2
src/include/PMsrHandler.h
View File

@ -107,7 +107,7 @@ class PMsrHandler
virtual void GetGroupingString(Int_t runNo, TString detector, TString &groupingStr);
virtual Bool_t EstimateN0();
virtual Double_t GetAlphaEstimateN0();
//as virtual Double_t GetAlphaEstimateN0();
private:
Bool_t fFourierOnly; ///< flag indicating if Fourier transform only is wished. If yes, some part of the msr-file blocks are not needed.

6
src/include/PMusrCanvas.h
View File

@ -206,12 +206,12 @@ class PMusrCanvas : public TObject, public TQObject
PMusrCanvas();
PMusrCanvas(const Int_t number, const Char_t* title,
Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, const Bool_t batch,
const Bool_t fourier=false);
const Bool_t fourier=false, const Bool_t avg=false);
PMusrCanvas(const Int_t number, const Char_t* title,
Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh,
PMsrFourierStructure fourierDefault,
const PIntVector markerList, const PIntVector colorList, const Bool_t batch,
const Bool_t fourier=false);
const Bool_t fourier=false, const Bool_t avg=false);
virtual ~PMusrCanvas();
virtual Bool_t IsValid() { return fValid; }
@ -236,6 +236,7 @@ class PMusrCanvas : public TObject, public TQObject
private:
Bool_t fStartWithFourier; ///< flag if true, the Fourier transform will be presented bypassing the time domain representation
Bool_t fStartWithAvg; ///< flag if true, the averaged data/Fourier will be presented
Int_t fTimeout; ///< timeout after which the Done signal should be emited. If timeout <= 0, no timeout is taking place
Bool_t fScaleN0AndBkg; ///< true=N0 and background is scaled to (1/ns), otherwise (1/bin) for the single histogram case
Bool_t fBatchMode; ///< musrview in ROOT batch mode
@ -312,6 +313,7 @@ class PMusrCanvas : public TObject, public TQObject
virtual void CleanupFourierDifference();
virtual void CleanupAverage();
virtual void CalcPhaseOptReFT();
virtual Double_t CalculateDiff(const Double_t x, const Double_t y, TH1F *theo);
virtual Double_t CalculateDiff(const Double_t x, const Double_t y, TGraphErrors *theo);
virtual Int_t FindBin(const Double_t x, TGraphErrors *graph);

3
src/include/PRunListCollection.h
View File

@ -76,6 +76,9 @@ 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 Double_t GetSingleHistoMaximumLikelihoodExpected(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetSingleRunMaximumLikelihood(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const;
virtual UInt_t GetTotalNoOfBinsFitted() const;

1
src/include/PRunSingleHisto.h
View File

@ -45,6 +45,7 @@ class PRunSingleHisto : public PRunBase
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 Double_t CalcMaxLikelihoodExpected(const std::vector<Double_t>& par);
virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins();