[ROOT] Bug on TTree::UnbinnedFit?

From: Hajime Nanjyo (nanjyo@icepp.s.u-tokyo.ac.jp)
Date: Fri Sep 10 2004 - 17:04:48 MEST


Dear ROOTers,

I found a bug in the function TreeUnbinnedFitLikelihood in TTreePlayer.cxx,
which was related to a bug in TTree::UnbinnedFit.

In the function TreeUnbinnedFitLikelihood, event weights were treated as
follows.
prob = fitfunc->EvalPar(x,par) * weight[i]/sum;

I think it should be
prob = TMath::Power(fitfunc->EvalPar(x,par),weight[i]/sum).

So I propose to change the original code,

     prob = fitfunc->EvalPar(x,par) * weight[i]/sum;
     if(prob > 0) logL += TMath::Log(prob);
     else         logL += logEpsilon;

to

     prob = fitfunc->EvalPar(x,par);
     if(prob > 0) logL += TMath::Log(prob) * weight[i]/sum;
     else         logL += logEpsilon * weight[i]/sum;
.

I show a small sample, "test.C" to indicate the problem below.

It can be compiled and executed as follows.
g++ -o test `root-config --cflags` `root-config --libs` -lMinuit test.C
./test

In this code, 20000 events with weights are generated and fitted with two ways.
1. Fill a histogram with the weights and fit it.
2. TTree::UnbinnedFit with weights.

The results should be same but not the same currently.
I already checked that the fix of the TreeUnbinnedFitLikelihood
remedy the problem.

The ROOT under /afs/cern.ch/sw/root/v4.00.08/rh73_gcc32/root was used 
with gcc version 3.2 in Red Hat Linux release 7.3 (Valhalla)

Could you treat the problem please.

Best Regards,
Hajime


test.C
//////////////////////////////////////////////////////////////
#include <iostream>
#include <cstdio>
#include "TROOT.h"
#include "TApplication.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TH1.h"

int main(int argc,char** argv)
{
  TApplication theApp("App", &argc, argv);
  gApplication->Init();

  //---------------------------------
  // 1.) generate n0 events following 1+x^2+8/3*a0*x with weight w0.
  // 2.) generate n1 events following 1+x^2+8/3*a1*x with weight w1.
  // 3.) fill (n0+n1) events in TTree tr0 in TFile test.root.
  // 4.) Fill values to the TH1D his with weight.
  // 5.) Fit the histogram with [0]*(1+x^2+8/3*[1]*x)
  // 6.) Performe UnbinnedFit with p.d.f of 1+x^2+8/3*[0]*x
  //---------------------------------
  
  int    n0=10000;
  int    n1=10000;

  double a0=0.3;
  double a1=0.6;
  double w0=0.4;
  double w1=0.8;
  //---------------------------------

  
  TFile tf("test.root","RECREATE");
  TF1 f0("f0","1+x*x+8./3.*[0]*x",-1,1);
  f0.SetParameter(0,a0);
  f0.SetNpx(1000);
  TF1 f1("f1","1+x*x+8./3.*[0]*x",-1,1);
  f1.SetParameter(0,a1);
  f1.SetNpx(1000);

  Double_t value, weight;
  
  TTree tr0("tr0","tr0");
  tr0.Branch("value",&value,"value/D");
  tr0.Branch("weight",&weight,"weight/D");
  
  for(int i=0;i<n0;i++) {
    value=f0.GetRandom();
    weight=w0;
    tr0.Fill();
  }
  for(int i=0;i<n1;i++) {
      value=f1.GetRandom();
    weight=w1;
    tr0.Fill();
  }
  tr0.Write();
  tf.Close();
  
  //---------------------------------
  
  TFile tfi("test.root");
  TTree* tr;
  tr = (TTree*)tfi.Get("tr0");
  
  TF1 fun("fun","[0]*(1+x*x+8./3.*[1]*x)",-1,1);

  TH1D his("his","his",40,-1,1);
  tr->Draw("value>>his","weight");
  his.Fit("fun","","",-1,1);
  Double_t A0= fun.GetParameter(1);
  //his.Draw("E");
  //gPad->Modified();
  //gPad->Update();
  //std::cout << "Hit Return" << std::endl;
  //std::getchar();
  //---------------------------------
  
  TF1 func("func","1+x*x+8./3.*[0]*x",-1,1);
  tr->UnbinnedFit("func","value","weight","VEM");
  Double_t A1= func.GetParameter(0);
  
  tfi.Close();

  std::cout << A0 << "\t" << A1 << std::endl;
  return 0;
}
//////////////////////////////////////////////////////////////  



This archive was generated by hypermail 2b29 : Sun Jan 02 2005 - 05:50:09 MET