Re: [ROOT] Bug on TTree::UnbinnedFit?

From: Rene Brun (Rene.Brun@cern.ch)
Date: Fri Sep 10 2004 - 18:46:41 MEST


Hi Hajime,

Thanks very much for your correction when using the weights.
iIt clearly improves the fit.

It is now in CVS.

Rene Brun

Hajime Nanjyo wrote:
> 
> 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