#include #include #include "TMath.h" #include "TGraph.h" #include "TF1.h" #include "TRandom.h" #include "Math/RootFinder.h" #include "Math/Functor.h" // Find multiple root of a function // script must be run in compiled mode for instantiating correctly the GradFunctor1D object // .x FindRoot.C+ bool IsEqual(double a, double b, double eps) { return TMath::Abs(a-b) < eps * TMath::Max( TMath::Abs(a), TMath::Abs(b) ); } void Solve(ROOT::Math::RootFinder & rf, vector& roots, int nest, const ROOT::Math::IGradFunction & f, double s, double e) { double x0 = (s+e)/2; rf.SetFunction(f, x0); double tol = 1.E-6; int status = rf.Solve(100,tol,tol); double r = rf.Root(); // is status is ok and root is in the interval if ( status == 0 && r > s && r < e) { cout << "(" << nest << ") start=" << s << " end=" << e << " root=" << r << endl; roots.push_back(r); if ( !IsEqual(r,s,tol) && ! IsEqual(r,e,tol ) ) { Solve(rf, roots, nest+1, f, s, r); Solve(rf, roots, nest+1, f, r, e); } } else { std::cout << "(" << nest << ") start=" << s << " end=" << e; if (status != 0) std::cout << " Error solving for roots - status " << status << std::endl; else std::cout << " No root found in interval - r=" << r << std::endl; return; } } class MyFunctionObject { public: MyFunctionObject() { //fN = gRandom->Integer(10) + 1; fN = gRandom ->Poisson(8); fSigmas = new double[fN]; fMus = new double[fN]; fAmps = new double[fN]; cout << "We have " << fN << " pulses:" << endl; for(int i = 0; i < fN; ++i) { fSigmas[i] = TMath::Abs(gRandom->Gaus(3.12, 0.312)); fMus [i] = TMath::Abs(gRandom->Exp(3.7) + gRandom->Exp(8.8)); fAmps [i] = TMath::Abs(gRandom->Gaus(0.3, 0.14)); cout << fAmps[i] << " " << fMus[i] << " " << fSigmas[i] << endl; } } double operator()(double *x, double* par) { double r = 0; for(int i = 0; i < fN; ++i) { r += fAmps[i] * TMath::Gaus(x[0], fMus[i], fSigmas[i], kFALSE); } // Subtract in order to cross the y axis. double s = 0.1 * fN; // one third of the mean times the number of pulses. return r - s; } private: int fN; double* fSigmas; double* fMus; double* fAmps; }; void FindAllRoots() { double s = -50; double e = 50; cout << "start =" << s << endl; cout << "end =" << e << endl; // Create the function and wrap it MyFunctionObject* fobj = new MyFunctionObject(); TF1* f = new TF1("MyFunctionObject", fobj, s, e, 1, "MyFunctionObject"); f->Draw(); vector roots; ROOT::Math::RootFinder rf(ROOT::Math::RootFinder::kGSL_STEFFENSON); ROOT::Math::GradFunctor1D wf(*f); Solve(rf, roots, 0, wf, s, e); sort(roots.begin(), roots.end()); double* xx = new double[roots.size()]; double* yy = new double[roots.size()]; cout << "Found " << roots.size() << " roots:"; for (unsigned int i = 0; i < roots.size(); ++i) { xx [i] = roots[i]; yy [i] = 0; cout << xx[i] << " "; } cout << endl; TGraph* gr = new TGraph(roots.size(), xx, yy); gr->SetMarkerColor(kRed); gr->SetMarkerStyle(21); gr->Draw("P"); cout.flush(); }