From 4d9732faa195a02c5744b0150f82e22c3b9f6e64 Mon Sep 17 00:00:00 2001 From: nemu Date: Wed, 6 Feb 2008 09:41:49 +0000 Subject: [PATCH] implemented max log likelihood, and other things --- src/Makefile.msr2msr | 6 +-- src/ToDo.txt | 13 ++++--- src/classes/PFitter.cpp | 4 +- src/classes/PFitterFcn.cpp | 15 +++++++- src/classes/PMsrHandler.cpp | 2 +- src/classes/PRunListCollection.cpp | 20 +++++----- src/classes/PRunSingleHisto.cpp | 49 ++++++++++++++++++++++++- src/tests/rootMinuit2Test/minuit2test.C | 7 +++- 8 files changed, 89 insertions(+), 27 deletions(-) diff --git a/src/Makefile.msr2msr b/src/Makefile.msr2msr index 0dfa6b2f..4c5ab04e 100644 --- a/src/Makefile.msr2msr +++ b/src/Makefile.msr2msr @@ -60,9 +60,9 @@ LIBS = $(ROOTLIBS) -lXMLParser GLIBS = $(ROOTGLIBS) -lXMLParser # PSI libs -PSILIBS = -lTLemRunHeader -lPMusr +#PSILIBS = -lTLemRunHeader -lPMusr # Minuit2 lib -MNLIB = -L$(ROOTSYS)/lib -lMinuit2 +#MNLIB = -L$(ROOTSYS)/lib -lMinuit2 # Executable EXEC = msr2msr @@ -78,7 +78,7 @@ all: $(EXEC) $(EXEC): $(OBJS) @echo "---> Building $(EXEC) ..." /bin/rm -f $(SHLIB) - $(LD) $(OBJS) -o $(EXEC) $(GLIBS) $(PSILIBS) $(MNLIB) + $(LD) $(OBJS) -o $(EXEC) $(GLIBS) @echo "done" # clean up: remove all object file (and core files) diff --git a/src/ToDo.txt b/src/ToDo.txt index f121d901..162ce3ff 100644 --- a/src/ToDo.txt +++ b/src/ToDo.txt @@ -25,7 +25,7 @@ short term: * write a little standalone program which is converting the old msr-files to the new ones. **DONE** 08-01-04 -* implement max.likelihood fits +* implement max.likelihood fits **DONE** 08-02-06 (for Single Histo only. For the others not clear how to do it) * implement table based theory functions (LF stuff) @@ -33,14 +33,15 @@ short term: intermediate term: --------------------- -* start writing docu, i.e. transfering WKM doc from German -> English and +* start writing docu, i.e. transferring WKM doc from German -> English and describe new features with examples! +* introduce error numbers with corresponding docu/explanation in the latex-docu. * implement RRF stuff * implement NonMuSR stuff * implement access to user-function, i.e. functions not defined within musrfit, using the ROOT dictionary feature to look for them and use them. -* think about if it is worth to modifiy wkm in order to read +* think about if it is worth to modify wkm in order to read the new msr-files, so that wkmview can be used for an intermediate time. @@ -65,6 +66,8 @@ fixes: * Needed to change some part of the Minuit2 source code, otherwise I got a strange linker error message. In MnMinos.h the constructors (3 overloaded versions) are defined implicitly. When linking against libPMusr.so which is using libMinuit2Base.so, I got the error message from the linker that is cannot find the reference - to the second constructor. When transfering the implicit definition from the header file to the cxx, the + to the second constructor. When transferring the implicit definition from the header file to the cxx, the problem was gone. Since to fiddle in the Minuit2 source code directly is ugly, I should look for a way to - build libPMusr.so which doesn't have this problem. \ No newline at end of file + build libPMusr.so which doesn't have this problem. + + -> This problem has been resolved cleanly by using minuit2 delivered with root! \ No newline at end of file diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index b1ff98b1..44c3948e 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -178,7 +178,7 @@ bool PFitter::CheckCommands() continue; } else if (it->fLine.Contains("CHI_SQUARE")) { fUseChi2 = true; - } else if (it->fLine.Contains("MAX_LIKELYHOOD")) { + } else if (it->fLine.Contains("MAX_LIKELIHOOD")) { fUseChi2 = false; } else if (it->fLine.Contains("INTERACTIVE")) { fCmdList.push_back(PMN_INTERACTIVE); @@ -325,7 +325,7 @@ bool PFitter::ExecuteMinimize() ROOT::Minuit2::MnMinimize minimize((*fFitterFcn), fMnUserParams); // minimize - ROOT::Minuit2::FunctionMinimum min = minimize(); + ROOT::Minuit2::FunctionMinimum min = minimize(); if (!min.IsValid()) { cout << endl << "**WARNING**: PFitter::ExecuteMinimize(): Fit did not converge, sorry ..."; return false; diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index 5245bace..541c1e43 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -53,6 +53,19 @@ PFitterFcn::PFitterFcn(PRunListCollection *runList, bool useChi2) fUp = 0.5; fRunListCollection = runList; + + // check if max likelihood is used together with asymmetry/RRF/nonMusr data. + // if yes place a warning since this option is not implemented and a fall back + // to chi2 will be used. + if (!fUseChi2) { + if ((fRunListCollection->GetNoOfAsymmetry() > 0) || + (fRunListCollection->GetNoOfRRF() > 0) || + (fRunListCollection->GetNoOfNonMusr() > 0)) { + cout << endl << "**WARNING**: Maximum Log Likelihood Fit is only implemented for Single Histogram Fit"; + cout << endl << " Will fall back to Chi Square Fit."; + cout << endl; + } + } } //-------------------------------------------------------------------------- @@ -80,7 +93,7 @@ double PFitterFcn::operator()(const std::vector& par) const value += fRunListCollection->GetAsymmetryChisq(par); value += fRunListCollection->GetRRFChisq(par); value += fRunListCollection->GetNonMusrChisq(par); - } else { // max likelyhood + } else { // max likelihood value += fRunListCollection->GetSingleHistoMaximumLikelihood(par); value += fRunListCollection->GetAsymmetryMaximumLikelihood(par); value += fRunListCollection->GetRRFMaximumLikelihood(par); diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index fd758c8e..55f50069 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1467,7 +1467,7 @@ bool PMsrHandler::FilterFunMapNumber(TString str, const char *filter, int &no) //-------------------------------------------------------------------------- /** *

In the command block there shall be a new command which can be used - * to switch between chi2 and maximum_likelyhood!! + * to switch between chi2 and maximum_likelihood!! * * \param lines is a list of lines containing the fitparameter block */ diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 972fd54b..fdc99923 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -99,9 +99,6 @@ bool PRunListCollection::Add(int runNo) int fitType = (*fMsrInfo->GetMsrRunList())[runNo].fFitType; -PRunAsymmetry* asym; -bool status; - switch (fitType) { case PRUN_SINGLE_HISTO: fRunSingleHistoList.push_back(new PRunSingleHisto(fMsrInfo, fData, runNo)); @@ -110,8 +107,6 @@ bool status; break; case PRUN_ASYMMETRY: fRunAsymmetryList.push_back(new PRunAsymmetry(fMsrInfo, fData, runNo)); - asym = fRunAsymmetryList[fRunAsymmetryList.size()-1]; - status = asym->IsValid(); if (!fRunAsymmetryList[fRunAsymmetryList.size()-1]->IsValid()) success = false; break; @@ -217,14 +212,15 @@ double PRunListCollection::GetSingleHistoMaximumLikelihood(const std::vector + *

Since it is not clear yet how to handle asymmetry fits with max likelihood + * the chi square will be used! */ double PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector& par) { double mlh = 0.0; for (unsigned int i=0; iCalcMaxLikelihood(par); + mlh += fRunAsymmetryList[i]->CalcChiSquare(par); return mlh; } @@ -233,14 +229,15 @@ double PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector + *

Since it is not clear yet how to handle RRF fits with max likelihood + * the chi square will be used! */ double PRunListCollection::GetRRFMaximumLikelihood(const std::vector& par) { double mlh = 0.0; for (unsigned int i=0; iCalcMaxLikelihood(par); + mlh += fRunRRFList[i]->CalcChiSquare(par); return mlh; } @@ -249,14 +246,15 @@ double PRunListCollection::GetRRFMaximumLikelihood(const std::vector& pa // GetNonMusrMaximumLikelihood //-------------------------------------------------------------------------- /** - *

+ *

Since it is not clear yet how to handle non musr fits with max likelihood + * the chi square will be used! */ double PRunListCollection::GetNonMusrMaximumLikelihood(const std::vector& par) { double mlh = 0.0; for (unsigned int i=0; iCalcMaxLikelihood(par); + mlh += fRunNonMusrList[i]->CalcChiSquare(par); return mlh; } diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 7571e31c..3c288b75 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -166,9 +166,54 @@ double PRunSingleHisto::CalcChiSquare(const std::vector& par) */ double PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) { - cout << endl << "PRunSingleHisto::CalcMaxLikelihood(): not implemented yet ..." << endl; + double mllh = 0.0; // maximum log likelihood assuming poisson distribution for the single bin - return 1.0; + double N0; + + // check if norm is a parameter or a function + if (fRunInfo->fNormParamNo < MSR_PARAM_FUN_OFFSET) { // norm is a parameter + N0 = par[fRunInfo->fNormParamNo-1]; + } else { // norm is a function + // get function number + unsigned int funNo = fRunInfo->fNormParamNo-MSR_PARAM_FUN_OFFSET; + // evaluate function + N0 = fMsrInfo->EvalFunc(funNo,fRunInfo->fMap,par); + } + + // get tau + double tau; + if (fRunInfo->fLifetimeParamNo != -1) + tau = par[fRunInfo->fLifetimeParamNo-1]; + else + tau = PMUON_LIFETIME; + + // get background + double bkg = par[fRunInfo->fBkgFitParamNo-1]; + + // calculate functions + for (int i=0; iGetNoOfFuncs(); i++) { + int funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, fRunInfo->fMap, par); + } + + // calculate maximum log likelihood + double theo; + double data; + for (unsigned int i=0; i=fFitStartTime) && (fData.fTime[i]<=fFitStopTime)) { + // calculate theory for the given parameter set + theo = N0*TMath::Exp(-fData.fTime[i]/tau)*(1+fTheory->Func(fData.fTime[i], par, fFuncValues))+bkg; + // check if data value is not too small + if (fData.fValue[i] > 1.0e-9) + data = fData.fValue[i]; + else + data = 1.0e-9; + // add maximum log likelihood contribution of bin i + mllh -= data*TMath::Log(theo) - theo - TMath::LnGamma(data+1); + } + } + + return mllh; } //-------------------------------------------------------------------------- diff --git a/src/tests/rootMinuit2Test/minuit2test.C b/src/tests/rootMinuit2Test/minuit2test.C index bbeaa56e..1a6b4eb5 100644 --- a/src/tests/rootMinuit2Test/minuit2test.C +++ b/src/tests/rootMinuit2Test/minuit2test.C @@ -79,7 +79,7 @@ void minuit2test() gFitFcn->SetParNames("N0", "asym", "lambda", "B", "phase", "Bkg"); gFitFcn->SetParameter(0, 30.0); // N0 // gFitFcn->SetParLimits(0, 0.0, 1.0e6); - gFitFcn->SetParameter(1, 0.11); // asym + gFitFcn->SetParameter(1, 0.26); // asym // gFitFcn->SetParLimits(1, 0.0, 0.33); gFitFcn->SetParameter(2, 0.3); // lambda // gFitFcn->SetParLimits(2, 0.0, 100.0); @@ -100,6 +100,9 @@ void minuit2test() histo->Draw(); + gFitFcn->SetParameter(0, 150.0); // N0 + gFitFcn->SetParameter(5, 23.0); // Bkg + TVirtualFitter::SetDefaultFitter("Minuit2"); - histo->Fit("gFitFcn", "LE"); // L->likleyhood, E->minos + histo->Fit("gFitFcn", "L"); // L->likleyhood, E->minos }