testNdimFit.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 = 10; 
const std::string branchType = "x[10]/D";
const int NPoints = 100000;

// const int N = 50; 
// const std::string branchType = "x[50]/D";
// const int NPoints = 10000;



double truePar[2*N]; 
double iniPar[2*N]; 
const int nfit = 1;
const int strategy = 0;

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

   double invsig = 1./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 gausnormN(const double *x, const double *p) { 
   double f = 1; 
   for (int i = 0; i < N; ++i) 
      f *= gausnorm(x+i,p+2*i);

   return f; 
}

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
   
   delete tree; 
   tree = 0; 
   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 << "data size type and size  is " << typeid(*d).name() <<  "   " << d->Size() << std::endl;
   std::cout << "Fit data size =  " << d->Size() << " dimension = " << d->NDim() << std::endl;

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

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

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


   // create the function

   fitter.SetFunction(func); 
   // need to fix param 0 , normalization in the unbinned fits
   //fitter.Config().ParSettings(0).Fix();

   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;
   for (int i = 0; i < N; ++i) { 
      double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
      chi2 += d*d; 
   }
   double prob = ROOT::Math::chisquared_cdf_c(chi2,N);
   int iret =  (prob < 1.0E-6) ? -1 : 0;
   if (iret != 0) {
      std::cout <<"Found difference in fitted values - prob = " << prob << std::endl;
      if (!debug) fitter.Result().Print(std::cout);    
      for (int i = 0; i < N; ++i) { 
         double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
         std::cout << "par_" << i << " = " << fitter.Result().Value(i) << " true  = " << truePar[i] << " pull = " << d << std::endl;
      } 
      
   }

   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 );

#else
   for (int i = 0; i < nfit; ++i) { 
      func.SetParameters(iniPar);
      iret = DoFit<MinType>(fitobj,func, false);
      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; 
}


int testNdimFit() { 


   std::cout << "\n\n************************************************************\n"; 
   std::cout << "\t UNBINNED TREE (GAUSSIAN MULTI-DIM)  FIT\n";
   std::cout << "************************************************************\n"; 

   TTree * t1 = new TTree("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 < N; ++j) {          
         double mu = double(j)/10.; 
         double s  = 1.0 + double(j)/10.;  
         truePar[2*j] = mu; 
         truePar[2*j+1] = s;
      }

   
   //fill the tree
   TRandom3 r; 
   for (Int_t i=0;i<NPoints;i++) {
      for (int j = 0;  j < N; ++j) { 
         double mu = truePar[2*j]; 
         double s  = truePar[2*j+1];  
         x[j] = r.Gaus(mu,s);
      }

      ev = i;
      t1->Fill();
      
   }
   //t1.Draw("x"); // to select fit variable 


   for (int i = 0; i <N; ++i) {
      iniPar[2*i] = 0; 
      iniPar[2*i+1] = 1; 
   }

   // use simply TF1 wrapper 
   //ROOT::Math::WrappedMultiTF1 f2(*f1); 
   ROOT::Math::WrappedParamFunction<> f2(&gausnormN,N,2*N,iniPar); 

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

   return iret; 
}

int main() { 
  return testNdimFit();
}

---