From c768c27898f93e7271b0195d438342718278cf02 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 18 Dec 2016 10:36:02 +0100 Subject: [PATCH] for maxLH it is now possible to get the per-run-block maxLH, and the expected maxLH --- src/classes/PFitter.cpp | 104 +++++----- src/classes/PFitterFcn.cpp | 11 +- src/classes/PMsrHandler.cpp | 305 ++++++++++++++++------------- src/classes/PRunListCollection.cpp | 93 ++++++++- src/classes/PRunSingleHisto.cpp | 100 +++++++++- src/include/PMsrHandler.h | 2 +- src/include/PMusrCanvas.h | 6 +- src/include/PRunListCollection.h | 3 + src/include/PRunSingleHisto.h | 1 + 9 files changed, 433 insertions(+), 192 deletions(-) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index e366bda1..c921c250 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -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 param; - Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerHisto; - std::vector ndfPerHisto; + // calculate expected chisq + std::vector param; + Double_t totalExpectedChisq = 0.0; + std::vector expectedChisqPerHisto; + std::vector ndfPerHisto; - for (UInt_t i=0; iCalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + // CalcExpectedChiSquare handles both, chisq and mlh + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); - // calculate chisq per run - std::vector chisqPerHisto; - for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + // calculate chisq per run + std::vector chisqPerHisto; + for (UInt_t i=0; iGetMsrRunList()->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; iGetNoOfFittedBins(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; iGetNoOfFittedBins(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; iGetNoOfFittedBins(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; iGetNoOfFittedBins(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; diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index c6a8566a..07c4831a 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -110,15 +110,20 @@ void PFitterFcn::CalcExpectedChiSquare(const std::vector &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; iGetNoOfSingleHisto(); 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; iGetNoOfSingleHisto(); i++) { + value = fRunListCollection->GetSingleHistoMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i' + expectedChisqPerRun.push_back(value); + totalExpectedChisq += value; + } } } diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index d126ccd9..dedf2d29 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -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; jsize(); 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 0) { + for (UInt_t i=0; 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; iwriteExpectedChisq) 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; iwriteExpectedChisq) + 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 0) { + for (UInt_t i=0; 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; iwriteExpectedChisq) 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; iwriteExpectedChisq) + 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 0) { + for (UInt_t i=0; 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; iwriteExpectedChisq) 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; iwriteExpectedChisq) + 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 0) { + for (UInt_t i=0; 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; iwriteExpectedChisq) 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; iwriteExpectedChisq) + 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() /** *

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) diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index c9c5e56f..3cc36d01 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -355,7 +355,7 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vectorGetMsrRunList()->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::vectorCalculates expected mlh of the single histogram with run block index idx of a msr-file. + * + * return: + * - expected mlh of for a single histogram + * + * \param par fit parameter vector + * \param idx run block index + */ +Double_t PRunListCollection::GetSingleHistoMaximumLikelihoodExpected(const std::vector& 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; iGetMsrRunList()->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) +//-------------------------------------------------------------------------- +/** + *

Calculates mlh of a single run-block entry of the msr-file. + * + * return: + * - 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& 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; iGetMsrRunList()->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) //-------------------------------------------------------------------------- diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 8b8b3094..23738b70 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -393,6 +393,103 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) return 2.0*mllh; } +//-------------------------------------------------------------------------- +// CalcMaxLikelihoodExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected log maximum-likelihood. + * + * return: + * - log maximum-likelihood value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector& 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; iGetNoOfFuncs(); 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; iFunc(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 0) denom += xx*xx/dval; diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index 1dab817b..097ccdc1 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -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. diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index cdc6791a..d61392c2 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -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); diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 00b1112e..1c22ee51 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -76,6 +76,9 @@ class PRunListCollection virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; + virtual Double_t GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleRunMaximumLikelihood(const std::vector& par, const UInt_t idx) const; + virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const; virtual UInt_t GetTotalNoOfBinsFitted() const; diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 15423636..ba8b49d4 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -45,6 +45,7 @@ class PRunSingleHisto : public PRunBase virtual Double_t CalcChiSquare(const std::vector& par); virtual Double_t CalcChiSquareExpected(const std::vector& par); virtual Double_t CalcMaxLikelihood(const std::vector& par); + virtual Double_t CalcMaxLikelihoodExpected(const std::vector& par); virtual void CalcTheory(); virtual UInt_t GetNoOfFitBins();