implemented max log likelihood, and other things

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)

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.
+  
+  -> This problem has been resolved cleanly by using minuit2 delivered with root!

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);

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<double>& 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);

src/classes/PMsrHandler.cpp

@@ -1467,7 +1467,7 @@ bool PMsrHandler::FilterFunMapNumber(TString str, const char *filter, int &no)
 //--------------------------------------------------------------------------
 /**
  * <p> 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
  */

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<dou
 // GetAsymmetryMaximumLikelihood
 //--------------------------------------------------------------------------
 /**
- * <p>
+ * <p> 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<double>& par)
 {
   double mlh = 0.0;
 
   for (unsigned int i=0; i<fRunAsymmetryList.size(); i++)
-    mlh += fRunAsymmetryList[i]->CalcMaxLikelihood(par);
+    mlh += fRunAsymmetryList[i]->CalcChiSquare(par);
 
   return mlh;
 }
@@ -233,14 +229,15 @@ double PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector<doubl
 // GetRRFMaximumLikelihood
 //--------------------------------------------------------------------------
 /**
- * <p>
+ * <p> 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<double>& par)
 {
   double mlh = 0.0;
 
   for (unsigned int i=0; i<fRunRRFList.size(); i++)
-    mlh += fRunRRFList[i]->CalcMaxLikelihood(par);
+    mlh += fRunRRFList[i]->CalcChiSquare(par);
 
   return mlh;
 }
@@ -249,14 +246,15 @@ double PRunListCollection::GetRRFMaximumLikelihood(const std::vector<double>& pa
 // GetNonMusrMaximumLikelihood
 //--------------------------------------------------------------------------
 /**
- * <p>
+ * <p> 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<double>& par)
 {
   double mlh = 0.0;
 
   for (unsigned int i=0; i<fRunNonMusrList.size(); i++)
-    mlh += fRunNonMusrList[i]->CalcMaxLikelihood(par);
+    mlh += fRunNonMusrList[i]->CalcChiSquare(par);
 
   return mlh;
 }

src/classes/PRunSingleHisto.cpp

@@ -166,9 +166,54 @@ double PRunSingleHisto::CalcChiSquare(const std::vector<double>& par)
  */
 double PRunSingleHisto::CalcMaxLikelihood(const std::vector<double>& 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; i<fMsrInfo->GetNoOfFuncs(); 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<fData.fValue.size(); i++) {
+    if ((fData.fTime[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;
 }
 
 //--------------------------------------------------------------------------

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 
 }