Dear Rene, I think "logL += TMath::Log(prob) * weight[i]/sum;" should be "logL += TMath::Log(prob) * weight[i];". Don't divide "weight[i]" by "sum". I'm sorry . This makes difference in the error evaluation of the fit. I think the fix below is right... ////////////////////////////////////////////////////////// prob = fitfunc->EvalPar(x,par); if(prob > 0) logL += TMath::Log(prob) * weight[i]; else logL += logEpsilon * weight[i]; ////////////////////////////////////////////////////////// Best Regards, Hajime > 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