#include "TH1.h" #include "TF1.h" #include "TF2.h" #include "TMath.h" #include "TSystem.h" #include "TRandom3.h" #include "TTree.h" #include "TROOT.h" #include "Fit/BinData.h" #include "Fit/UnBinData.h" //#include "Fit/BinPoint.h" #include "Fit/Fitter.h" #include "HFitInterface.h" #include "Math/IParamFunction.h" #include "Math/WrappedTF1.h" #include "Math/WrappedMultiTF1.h" #include "Math/WrappedParamFunction.h" #include "Math/MultiDimParamFunctionAdapter.h" #include "TGraphErrors.h" #include "TStyle.h" #include "TSeqCollection.h" #include "Math/Polynomial.h" #include "Math/DistFunc.h" #include #include #include "TStopwatch.h" #include "TVirtualFitter.h" // #include "TFitterFumili.h" // #include "TFumili.h" #include "GaussFunction.h" // #include "RooDataHist.h" // #include "RooDataSet.h" // #include "RooRealVar.h" // #include "RooGaussian.h" // #include "RooMinuit.h" // #include "RooChi2Var.h" // #include "RooGlobalFunc.h" // #include "RooFitResult.h" // #include "RooProdPdf.h" #include #include "MinimizerTypes.h" #ifdef USE_VC #include "Vc/Vc" //#include "Vc/Allocator" //Vc_DECLARE_ALLOCATOR(Vc::double_v) #endif //#define USE_AOS #ifdef USE_VDT #include "vdtMath.h" #endif //#define DEBUG int nfit; const int N = 20; double iniPar[2*N]; int ndimObs; int ndimPars; void printData(const ROOT::Fit::UnBinData & data) { for (unsigned int i = 0; i < data.Size(); ++i) { std::cout << data.Coords(i)[0] << "\t"; } std::cout << "\ndata size is " << data.Size() << std::endl; } void printResult(int iret) { std::cout << "\n************************************************************\n"; std::cout << "Test\t\t\t\t"; if (iret == 0) std::cout << "OK"; else std::cout << "FAILED"; std::cout << "\n************************************************************\n"; } bool USE_BRANCH = false; ROOT::Fit::UnBinData * FillUnBinData(TTree * tree, bool copyData = true, unsigned int dim = 1 ) { // fill the unbin data set from a TTree ROOT::Fit::UnBinData * d = 0; // for the large tree if (std::string(tree->GetName()) == "t2") { d = new ROOT::Fit::UnBinData(); // large tree unsigned int n = tree->GetEntries(); #ifdef DEBUG std::cout << "number of unbin data is " << n << " of dim " << N << std::endl; #endif d->Initialize(n,N); TBranch * bx = tree->GetBranch("x"); double vx[N]; bx->SetAddress(vx); std::vector m(N); for (int unsigned i = 0; i < n; ++i) { bx->GetEntry(i); d->Add(vx); for (int j = 0; j < N; ++j) m[j] += vx[j]; } #ifdef DEBUG std::cout << "average values of means :\n"; for (int j = 0; j < N; ++j) std::cout << m[j]/n << " "; std::cout << "\n"; #endif return d; } if (USE_BRANCH) { d = new ROOT::Fit::UnBinData(); unsigned int n = tree->GetEntries(); //std::cout << "number of unbin data is " << n << std::endl; if (dim == 2) { d->Initialize(n,2); TBranch * bx = tree->GetBranch("x"); TBranch * by = tree->GetBranch("y"); double v[2]; bx->SetAddress(&v[0]); by->SetAddress(&v[1]); for (int unsigned i = 0; i < n; ++i) { bx->GetEntry(i); by->GetEntry(i); d->Add(v); } } else if (dim == 1) { d->Initialize(n,1); TBranch * bx = tree->GetBranch("x"); double v[1]; bx->SetAddress(&v[0]); for (int unsigned i = 0; i < n; ++i) { bx->GetEntry(i); d->Add(v); } } return d; //printData(d); } else { tree->SetEstimate(tree->GetEntries()); // use TTREE::Draw if (dim == 2) { tree->Draw("x:y",0,"goff"); // goff is used to turn off the graphics double * x = tree->GetV1(); double * y = tree->GetV2(); if (x == 0 || y == 0) { USE_BRANCH= true; return FillUnBinData(tree, true, dim); } // use array pre-allocated in tree->Draw . This is faster //assert(x != 0); unsigned int n = tree->GetSelectedRows(); if (copyData) { d = new ROOT::Fit::UnBinData(n,2); double vx[2]; for (int unsigned i = 0; i < n; ++i) { vx[0] = x[i]; vx[1] = y[i]; d->Add(vx); } } else // use data pointers directly d = new ROOT::Fit::UnBinData(n,x,y); } else if ( dim == 1) { tree->Draw("x",0,"goff"); // goff is used to turn off the graphics double * x = tree->GetV1(); if (x == 0) { USE_BRANCH= true; return FillUnBinData(tree, true, dim); } unsigned int n = tree->GetSelectedRows(); if (copyData) { d = new ROOT::Fit::UnBinData(n,1); for (int unsigned i = 0; i < n; ++i) { d->Add(x[i]); } } else d = new ROOT::Fit::UnBinData(n,x); } return d; } //std::copy(x,x+n, d.begin() ); return 0; } // print the data template void printData(const T & data) { for (typename T::const_iterator itr = data.begin(); itr != data.end(); ++itr) { std::cout << itr->Coords()[0] << " " << itr->Value() << " " << itr->Error() << std::endl; } std::cout << "\ndata size is " << data.Size() << std::endl; } // new likelihood function for unbinned data using namespace ROOT::Fit; template struct VecNLL { VecNLL( ROOT::Fit::UnBinData & data, F & f) : fData(&data), fFunc(f) { } //GetData(); } // void GetData() { // const UnBinData & data = *fData; // //const ROOT::Math::IParamMultiFunction & func = *fFunc; //fNLL.ModelFunction(); // unsigned int n = data.Size(); // fX = std::vector(n); // fY = std::vector(n); // for (int i = 0; iGetNpar()); pol.SetParameters(func->GetParameters() ); ROOT::Math::WrappedParamFunction f(pol,1,func->GetParameters(),func->GetParameters()+func->GetNpar() ); #endif GaussFunction f; f.SetParameters(func.Parameters()); fitter.SetFunction(f); } bool ret = fitter.Fit(d); if (!ret) { std::cout << " Fit Failed " << std::endl; return -1; } if (debug) fitter.Result().Print(std::cout); return 0; } // unbin fit template int DoUnBinFit(T * tree, Func & func, bool debug = false, bool copyData = false ) { ROOT::Fit::UnBinData * d = FillUnBinData(tree, copyData, func.NDim() ); // need to have done Tree->Draw() before fit //FillUnBinData(d,tree); //std::cout << "data size type and size is " << typeid(*d).name() << " " << d->Size() << std::endl; if (debug) { if (copyData) std::cout << "\tcopy data in FitData\n"; else std::cout << "\tre-use original data \n"; } //printData(d); // create the fitter //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl; ROOT::Fit::Fitter fitter; fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str()); if (debug) fitter.Config().MinimizerOptions().SetPrintLevel(3); // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt fitter.Config().MinimizerOptions().SetTolerance(1); // create the function // need to fix param 0 , normalization in the unbinned fits //fitter.Config().ParSettings(0).Fix(); fitter.SetFunction(func); bool ret = fitter.Fit(*d); if (!ret) { std::cout << " Fit Failed " << std::endl; return -1; } if (debug) fitter.Result().Print(std::cout); delete d; return 0; } template int DoUnBinFitVec(T * tree, VFunc & func, int ndim, int npar, const double * p0, bool debug = false, bool copyData = false ) { ROOT::Fit::UnBinData * d = FillUnBinData(tree, copyData, ndim ); // need to have done Tree->Draw() before fit //FillUnBinData(d,tree); //std::cout << "data size type and size is " << typeid(*d).name() << " " << d->Size() << std::endl; if (debug) { if (copyData) std::cout << "\tcopy data in FitData\n"; else std::cout << "\tre-use original data \n"; } //printData(d); // create the fitter //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl; ROOT::Fit::Fitter fitter; fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str()); if (debug) fitter.Config().MinimizerOptions().SetPrintLevel(3); // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt fitter.Config().MinimizerOptions().SetTolerance(1); // create the function // need to fix param 0 , normalization in the unbinned fits //fitter.Config().ParSettings(0).Fix(); VecNLL fcn(*d, func); bool ret = fitter.FitFCN(npar, fcn, p0); if (!ret) { std::cout << " Fit Failed " << std::endl; return -1; } if (debug) fitter.Result().Print(std::cout); std::cout << "use vec nll : nll = " << fitter.Result().MinFcnValue() << std::endl; delete d; return 0; } template int DoFit(TTree * tree, Func & func, bool debug = false, bool copyData = false ) { return DoUnBinFit(tree, func, debug, copyData); } template int DoFit(TH1 * h1, Func & func, bool debug = false, bool copyData = false ) { return DoBinFit(h1, func, debug, copyData); } template int DoFit(TGraph * gr, Func & func, bool debug = false, bool copyData = false ) { return DoBinFit(gr, func, debug, copyData); } template int DoFitVec(TTree * tree, F & func, int n1, int n2, const double * p, bool debug = false, bool copyData = false ) { return DoUnBinFitVec(tree, func, n1,n2,p,debug, copyData); } template int DoFitVec(TH1 * h1, F & func, int, int , const double *, bool debug = false, bool copyData = false ) { return DoBinFit(h1, func, debug, copyData); } template int DoFitVec(TGraph * gr, F & func, int, int , const double *, bool debug = false, bool copyData = false ) { return DoBinFit(gr, func, debug, copyData); } template int FitUsingNewFitter(FitObj * fitobj, FuncObj func, bool useGrad=false) { std::cout << "\n************************************************************\n"; std::cout << "\tFit using new Fit::Fitter " << typeid(*fitobj).name() << std::endl; std::cout << "\tMinimizer is " << MinType::name() << " " << MinType::name2() << " func dim = " << ndimObs << std::endl; int iret = 0; TStopwatch w; w.Start(); #define USE_VECNLL #ifndef USE_VECNLL #ifdef DEBUG // std::cout << "initial Parameters " << iniPar << " " << *iniPar << " " << *(iniPar+1) << std::endl; func.SetParameters(iniPar); iret |= DoFit(fitobj,func,true, useGrad); if (iret != 0) { std::cout << "Fit failed " << std::endl; } #else for (int i = 0; i < nfit; ++i) { func.SetParameters(iniPar); iret = DoFit(fitobj,func, false, useGrad); if (iret != 0) { std::cout << "Fit failed " << std::endl; break; } } #endif #else // use vectorized function for (int i = 0; i < nfit; ++i) { iret = DoFitVec(fitobj,func, ndimObs, ndimPars, iniPar, false, useGrad); if (iret != 0) { std::cout << "Fit failed " << std::endl; break; } } #endif w.Stop(); std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl; std::cout << "\n************************************************************\n"; return iret; } double poly2(const double *x, const double *p) { return p[0] + (p[1]+p[2]*x[0] ) * x[0]; } int testPolyFit() { int iret = 0; std::cout << "\n\n************************************************************\n"; std::cout << "\t POLYNOMIAL FIT\n"; std::cout << "************************************************************\n"; std::string fname("pol2"); //TF1 * func = (TF1*)gROOT->GetFunction(fname.c_str()); TF1 * f1 = new TF1("pol2",fname.c_str(),-5,5.); f1->SetParameter(0,1); f1->SetParameter(1,0.0); f1->SetParameter(2,1.0); // fill an histogram TH1D * h1 = new TH1D("h1","h1",20,-5.,5.); // h1->FillRandom(fname.c_str(),100); for (int i = 0; i <1000; ++i) h1->Fill( f1->GetRandom() ); //h1->Print(); //h1->Draw(); iniPar[0] = 2.; iniPar[1] = 2.; iniPar[2] = 2.; // dummy for testing //iret |= FitUsingNewFitter(h1,f1); // use simply TF1 wrapper //ROOT::Math::WrappedMultiTF1 f2(*f1); ROOT::Math::WrappedParamFunction<> f2(&poly2,1,iniPar,iniPar+3); // if Minuit2 is later than TMinuit on Interl is much slower , why ?? iret |= FitUsingNewFitter(h1,f2); iret |= FitUsingNewFitter(h1,f2); // test with linear fitter // for this test need to pass a multi-dim function ROOT::Math::WrappedTF1 wf(*f1); ROOT::Math::MultiDimParamGradFunctionAdapter lfunc(wf); iret |= FitUsingNewFitter(h1,lfunc,true); // test with a graph gStyle->SetErrorX(0.); // to seto zero error on X TGraphErrors * gr = new TGraphErrors(h1); iret |= FitUsingNewFitter(gr,f2); std::cout << "\n-----> test now TGraphErrors with errors in X coordinates\n\n"; // try with error in X gStyle->SetErrorX(0.5); // to set zero error on X TGraphErrors * gr2 = new TGraphErrors(h1); iret |= FitUsingNewFitter(gr2,f2); printResult(iret); return iret; } template struct GausFunctions { typedef T value_type; static T gaussian(const T *x, const double *p) { //return p[0]*TMath::Gaus(x[0],p[1],p[2]); T tmp = (x[0]-p[1])/p[2]; return p[0] * std::exp(-tmp*tmp/2); } static T gausnorm(const T *x, const double *p) { //return p[0]*TMath::Gaus(x[0],p[1],p[2]); T invsig = 1./p[1]; T tmp = (x[0]-p[0]) * invsig; const T sqrt_2pi = 1./std::sqrt(2.* 3.14159 ); return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig; } static T gausnorm2D(const T *x, const double *p) { //return p[0]*TMath::Gaus(x[0],p[1],p[2]); return gausnorm(x,p)*gausnorm(x+1,p+2); } static T gausnormN(const T *x, const double *p) { //return p[0]*TMath::Gaus(x[0],p[1],p[2]); T f = 1.0; for (int i = 0; i < N; ++i) f *= gausnorm(x+i,p+2*i); return f; } static void gausnorm_v(unsigned int n, unsigned int stride, const double *x, const double *p, double * res) { double invsig = 1./p[1]; std::vector arg(n); for (unsigned int i = 0; i< n; ++i){ double tmp = (x[i*stride]-p[0]) * invsig; arg[i] = -0.5 *tmp*tmp; } const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 ); #ifdef USE_VDT vdt::fast_expv(n,&arg[0],res); for (unsigned int i = 0; i< n; ++i){ #else for (unsigned int i = 0; i< n; ++i){ res[i] = std::exp(arg[i]); #endif res[i] *= sqrt_2pi * invsig; } } static void gausnorm2D_v(unsigned int n, unsigned int stride,const double *x, const double *p, double * res) { std::vector tmp(n); gausnorm_v(n,2,x,p,&tmp[0]); gausnorm_v(n,2,x+1,p+2,res); for (unsigned int i = 0; i< n; ++i){ res[i] *= tmp[i]; } } static void gausnormN_v(unsigned int n, unsigned int stride,const double *x, const double *p, double * res) { std::vector tmp(n); gausnorm_v(n,stride,x,p,res); for (int j = 1; j < stride; ++j) { gausnorm_v(n,stride,x+j,p+2*j,&tmp[0]); for (unsigned int i = 0; i< n; ++i){ res[i] *= tmp[i]; } } } }; double gaussian (const double *x, const double *p) { return GausFunctions::gaussian(x,p); } double gausnorm (const double *x, const double *p) { return GausFunctions::gausnorm(x,p); } double gausnorm2D (const double *x, const double *p) { return GausFunctions::gausnorm2D(x,p); } double gausnormN (const double *x, const double *p) { return GausFunctions::gausnormN(x,p); } int testGausFit() { int iret = 0; std::cout << "\n\n************************************************************\n"; std::cout << "\t GAUSSIAN FIT\n"; std::cout << "************************************************************\n"; //std::string fname = std::string("gaus"); //TF1 * func = (TF1*)gROOT->GetFunction(fname.c_str()); //TF1 * f1 = new TF1("gaus",fname.c_str(),-5,5.); //TF1 * f1 = new TF1("gaussian",gaussian,-5,5.,3); //f2->SetParameters(0,1,1); // fill an histogram int nbin = 10000; TH1D * h2 = new TH1D("h2","h2",nbin,-5.,5.); // h1->FillRandom(fname.c_str(),100); for (int i = 0; i < 10000000; ++i) h2->Fill( gRandom->Gaus(0,10) ); iniPar[0] = 100.; iniPar[1] = 2.; iniPar[2] = 2.; // use simply TF1 wrapper //ROOT::Math::WrappedMultiTF1 f2(*f1); ROOT::Math::WrappedParamFunction<> f2(&gaussian,1,iniPar,iniPar+3); iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); // iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); // test also fitting a TGraphErrors with histogram data gStyle->SetErrorX(0.); // to seto zero error on X TGraphErrors * gr = new TGraphErrors(h2); iret |= FitUsingNewFitter(gr,f2); // try with error in X gStyle->SetErrorX(0.5); // to seto zero error on X TGraphErrors * gr2 = new TGraphErrors(h2); iret |= FitUsingNewFitter(gr2,f2); //#ifdef LATER // test using grad function std::cout << "\n\nTest Using pre-calculated gradients\n\n"; bool useGrad=true; iret |= FitUsingNewFitter(h2,f2,useGrad); iret |= FitUsingNewFitter(h2,f2,useGrad); iret |= FitUsingNewFitter(h2,f2,useGrad); iret |= FitUsingNewFitter(h2,f2,useGrad); iret |= FitUsingNewFitter(h2,f2,useGrad); iret |= FitUsingNewFitter(h2,f2,useGrad); // test LS algorithm std::cout << "\n\nTest Least Square algorithms\n\n"; iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); iret |= FitUsingNewFitter(h2,f2); // iret |= FitUsingTFit(h2,f1); // iret |= FitUsingTFit(h2,f1); //#endif //iret |= FitUsingRooFit(h2,f1); printResult(iret); return iret; } int testTreeFit() { std::cout << "\n\n************************************************************\n"; std::cout << "\t UNBINNED TREE (GAUSSIAN) FIT\n"; std::cout << "************************************************************\n"; TTree t1("t1","a simple Tree with simple variables"); double x, y; Int_t ev; t1.Branch("x",&x,"x/D"); t1.Branch("y",&y,"y/D"); // t1.Branch("pz",&pz,"pz/F"); // t1.Branch("random",&random,"random/D"); t1.Branch("ev",&ev,"ev/I"); //fill the tree int nrows = 10000; #ifdef TREE_FIT2D nrows = 10000; #endif for (Int_t i=0;iRannor(x,y); x *= 2; x += 1.; y *= 3; y -= 2; ev = i; t1.Fill(); } //t1.Draw("x"); // to select fit variable //TF1 * f1 = new TF1("gausnorm", gausnorm, -10,10, 2); //TF2 * f2 = new TF2("gausnorm2D", gausnorm2D, -10,10, -10,10, 4); ROOT::Math::WrappedParamFunction<> wf1(&gausnorm,1,iniPar,iniPar+2); ROOT::Math::WrappedParamFunction<> wf2(&gausnorm2D,2,iniPar,iniPar+4); iniPar[0] = 0; iniPar[1] = 1; iniPar[2] = 0; iniPar[3] = 1; // use simply TF1 wrapper //ROOT::Math::WrappedMultiTF1 f2(*f1); int iret = 0; // fit 1D first // iret |= FitUsingNewFitter(&t1,wf1,false); // not copying the data // iret |= FitUsingNewFitter(&t1,wf1,false); // not copying the data ndimObs = wf1.NDim(); ndimPars = wf1.NPar(); #ifdef USE_AOS iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnorm_v,true); // copying the data #else #ifndef USE_VC iret |= FitUsingNewFitter(&t1,wf1,true); // copying the data iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnorm,true); // copying the data #else iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnorm,true); // copying the data #endif // iret |= FitUsingNewFitter(&t1,wf1,true); // copying the data #endif // fit 2D ndimObs = wf2.NDim(); ndimPars = wf2.NPar(); #ifndef USE_AOS #ifndef USE_VC iret |= FitUsingNewFitter(&t1,wf2, true); iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnorm2D, true); #else iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnorm2D, true); #endif #else iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnorm2D_v,true); // copying the data #endif //iret |= FitUsingNewFitter(&t1,wf2, false); printResult(iret); return iret; } int testLargeTreeFit(int nevt = 1000) { std::cout << "\n\n************************************************************\n"; std::cout << "\t UNBINNED TREE (GAUSSIAN MULTI-DIM) FIT\n"; std::cout << "************************************************************\n"; TTree t1("t2","a large Tree with simple variables"); double x[N]; Int_t ev; t1.Branch("x",x,"x[20]/D"); t1.Branch("ev",&ev,"ev/I"); //fill the tree TRandom3 r; for (Int_t i=0;i f2(&gausnormN,N,2*N,iniPar); ndimObs = f2.NDim(); ndimPars = f2.NPar(); int iret = 0; #ifndef USE_VC #ifndef USE_AOS iret |= FitUsingNewFitter(&t1,f2); iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnormN,true); #else iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnormN_v,true); #endif // iret |= FitUsingNewFitter(&t1,f2); #else iret |= FitUsingNewFitter(&t1,&GausFunctions::gausnormN,true); #endif printResult(iret); return iret; } int testFitPerf() { int iret = 0; #ifdef DEBUG nfit = 1; #else nfit = 10; #endif iret |= testTreeFit(); nfit = 1; iret |= testLargeTreeFit(2000); #ifdef LATER #ifndef DEBUG nfit = 10; #endif iret |= testGausFit(); #ifndef DEBUG nfit = 1000; #endif iret |= testPolyFit(); #endif //return iret; if (iret != 0) std::cerr << "testFitPerf :\t FAILED " << std::endl; else std::cerr << "testFitPerf :\t OK " << std::endl; return iret; } int main() { return testFitPerf(); }