testUnbinGausFit.cxx

Go to the documentation of this file.

#include "TMath.h"
#include "TSystem.h"
#include "TRandom3.h"
#include "TTree.h"
#include "TROOT.h"

#include "Fit/UnBinData.h"
#include "Fit/Fitter.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 "Math/DistFunc.h"

#include <string>
#include <iostream>

#include "TStopwatch.h"

#define DEBUG

// test fit with many dimension

const int N = 1; // 1d fit
const int NGaus = 3;
const int NPar = 8; // sum of 3 gaussians
const std::string branchType = "x[1]/D";

// for 8 core testing use 1M points
const int NPoints = 100000;
double truePar[NPar];
double iniPar[NPar];
//const int nfit = 1;
const int strategy = 0;

double gausnorm(const double *x, const double *p) {

   double invsig = 1./std::abs(p[1]);
   double tmp = (x[0]-p[0]) * invsig;
   const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 );
   return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig;
}

double gausSum(const double *x, const double *p) {

   double f = gausnorm(x,p+2) +
      p[0] *  gausnorm(x,p+4) +
      p[1]  * gausnorm(x,p+6);

   double norm = 1. + p[0] + p[1];
   return f/norm;
}

struct MINUIT2 {
   static std::string name() { return "Minuit2"; }
   static std::string name2() { return ""; }
};

// fill fit data structure
ROOT::Fit::UnBinData * FillUnBinData(TTree * tree ) {

   // fill the unbin data set from a TTree
   ROOT::Fit::UnBinData * d = 0;

   // for the large tree access directly the branches
   d = new ROOT::Fit::UnBinData();

   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<double>  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;
}


// unbin fit

typedef ROOT::Math::IParamMultiFunction Func;
template <class MinType, class T>
int DoUnBinFit(T * tree, Func & func, bool debug = false ) {

   ROOT::Fit::UnBinData * d  = FillUnBinData(tree );
   // need to have done Tree->Draw() before fit
   //FillUnBinData(d,tree);

   std::cout << "Filled the fit data " << std::endl;
   //printData(d);

#ifdef DEBUG
   std::cout << "data size type and size  is " << typeid(*d).name() <<  "   " << d->Size() << std::endl;
#endif



   // 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);
   else
      fitter.Config().MinimizerOptions().SetPrintLevel(1);


   // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt
   fitter.Config().MinimizerOptions().SetTolerance(0.01);

   // set strategy (0 to avoid MnHesse
   fitter.Config().MinimizerOptions().SetStrategy(strategy);


   // create the function

   fitter.SetFunction(func);
   // need to set limits to constant term
   fitter.Config().ParSettings(0).SetLowerLimit(0.);
   fitter.Config().ParSettings(1).SetLowerLimit(0.);

   if (debug)
     std::cout << "do fitting... " << std::endl;

   bool ret = fitter.Fit(*d);
   if (!ret) {
      std::cout << " Fit Failed " << std::endl;
      return -1;
   }
   if (debug)
      fitter.Result().Print(std::cout);

   // check fit result
   double chi2 = 0;
   //if (fitter.Result().Value(0) <  0.5 ) {
   for (int i = 0; i < NPar; ++i) {
      double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
      chi2 += d*d;
   }
//}
//    else {
//      double truePar2[NPar];
//      truePar2[0] = 1.-truePar[0];
//      truePar2[1] = truePar[3];
//      truePar2[2] = truePar[4];
//      truePar2[3] = truePar[1];
//      truePar2[4] = truePar[2];
//      for (int i = 0; i < N; ++i) {
//        double d = ( truePar2[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
//        chi2 += d*d;
//      }
//    }
   double prob = ROOT::Math::chisquared_cdf_c(chi2,NPar);
   int iret =  (prob < 1.0E-6) ? -1 : 0;
   if (iret != 0) {
      std::cout <<"Found difference in fitted values - chi2 = " << chi2
                << " prob = " << prob << std::endl;
      fitter.Result().Print(std::cout);
   }

   delete d;

   return iret;

}


template <class MinType>
int DoFit(TTree * tree, Func & func, bool debug = false ) {
   return DoUnBinFit<MinType, TTree>(tree, func, debug);
}
// template <class MinType>
// int DoFit(TH1 * h1, Func & func, bool debug = false, bool copyData = false ) {
//    return DoBinFit<MinType, TH1>(h1, func, debug, copyData);
// }
// template <class MinType>
// int DoFit(TGraph * gr, Func & func, bool debug = false, bool copyData = false ) {
//    return DoBinFit<MinType, TGraph>(gr, func, debug, copyData);
// }

template <class MinType, class FitObj>
int FitUsingNewFitter(FitObj * fitobj, Func & func ) {

   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 = " << func.NDim() << std::endl;

   int iret = 0;
   TStopwatch w; w.Start();

#ifdef DEBUG
   std::cout << "initial Parameters " << iniPar << "  " << *iniPar << "   " <<  *(iniPar+1) << std::endl;
   func.SetParameters(iniPar);
   iret |= DoFit<MinType>(fitobj,func,true );
   if (iret != 0) {
     std::cout << "Test  failed " << std::endl;
   }

#else
   for (int i = 0; i < nfit; ++i) {
      func.SetParameters(iniPar);
      iret = DoFit<MinType>(fitobj,func, false);
      if (iret != 0) {
         std::cout << "Test failed " << std::endl;
         break;
      }
   }
#endif
   w.Stop();
   std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
   std::cout << "\n************************************************************\n";

   return iret;
}


int testNdimFit() {


   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 gaussian variables");
   double  x[N];
   Int_t ev;
   t1.Branch("x",x,branchType.c_str());
   t1.Branch("ev",&ev,"ev/I");

   // generate the true parameters
//       for (int j = 0;  j < NGaus; ++j) {
//          double a = j+1;
//          double mu = double(j)/NGaus;
//          double s  = 1.0 + double(j)/NGaus;
//          truePar[3*j] = a;
//          truePar[3*j+1] = mu;
//          truePar[3*j+2] = s;
//          tot += a;
//       }
   truePar[0] = 0.2; // % second  gaussian
   truePar[1] = 0.05;  // % third gaussian ampl
   truePar[2] = 0.;      // mean first gaussian
   truePar[3] = 0.5;    // s1
   truePar[4] = 0.;   // mean secon gauss
   truePar[5] = 1;
   truePar[6] = -3;   // mean third gaus
   truePar[7] = 10;



   //fill the tree
   TRandom3 r;
   double norm = (1+truePar[0] + truePar[1] );
   double a = 1./norm;
   double b = truePar[0]/ norm;
   double c = truePar[1]/ norm;
   assert(a+b+c == 1.);
   std::cout << " True amplitude gaussians " << a << "  " << b << "  " << c << std::endl;
   for (Int_t i=0;i<NPoints;i++) {
      for (int j = 0;  j < N; ++j) {
         if (r.Rndm() < a ) {
            double mu = truePar[2];
            double s  = truePar[3];
            x[j] = r.Gaus(mu,s);
         }
         else if (r.Rndm() < b ) {
            double mu = truePar[4];
            double s  = truePar[5];
            x[j] = r.Gaus(mu,s);
         }
         else {
            double mu = truePar[6];
            double s  = truePar[7];
            x[j] = r.Gaus(mu,s);
         }
      }

      ev = i;
      t1.Fill();

   }
   //t1.Draw("x"); // to select fit variable

   iniPar[0] = 0.5;
   iniPar[1] = 0.05;
   for (int i = 0; i <NGaus; ++i) {
      iniPar[2*i+2] = 0 ;
      iniPar[2*i+3] = 1. + 4*i;
      std::cout << "inipar " << i << " = " << iniPar[2*i+2] << "  " << iniPar[2*i+3] << std::endl;
   }

   // use simply TF1 wrapper
   //ROOT::Math::WrappedMultiTF1 f2(*f1);
   ROOT::Math::WrappedParamFunction<> f2(&gausSum,1,NPar,iniPar);

   int iret = 0;
   iret |= FitUsingNewFitter<MINUIT2>(&t1,f2);

   return iret;
}

int main() {
  return testNdimFit();
}