add new tutorial showing a combined fit to two histograms

git-svn-id: http://root.cern.ch/svn/root/trunk@37450 27541ba8-7e3a-0410-8455-c3a389f83636

tutorials/fit/NumericalMinimization.C

-// tutorial showing how to use the Minimizer class in ROOT
-//  and all the possible minimizer 
+// Example on how to use the new Minimizer class in ROOT
+//  Show usage with all the possible minimizers. 
 // Minimize the Rosenbrock function (a 2D -function)
 // This example is described also in 
 // http://root.cern.ch/drupal/content/numerical-minimization#multidim_minim
 // input : minimizer name + algorithm name
 // randomSeed: = <0 : fixed value: 0 random with seed 0; >0 random with given seed 
+//
+//Author: L. Moneta Dec 2010
 
 #include "Math/Minimizer.h"
 #include "Math/Factory.h"
@@ -22,19 +24,24 @@ double RosenBrock(const double *xx )
   return 100*tmp1*tmp1+tmp2*tmp2;
 }
  
-int NumericalMinimization(const char * minName = "Minuit2",const char *algoName = "" , int randomSeed = -1)
+int NumericalMinimization(const char * minName = "Minuit2",
+                          const char *algoName = "" , 
+                          int randomSeed = -1)
 {
-   // create minimizer giving a name and a name (optionally) for the specific algorithm
+   // create minimizer giving a name and a name (optionally) for the specific
+   // algorithm
    // possible choices are: 
    //     minName                  algoName
    // Minuit /Minuit2             Migrad, Simplex,Combined,Scan  (default is Migrad)
    //  Minuit2                     Fumili2
    //  Fumili
-   //  GSLMultiMin                ConjugateFR, ConjugatePR, BFGS, BFGS2, SteepestDescent
+   //  GSLMultiMin                ConjugateFR, ConjugatePR, BFGS, 
+   //                              BFGS2, SteepestDescent
    //  GSLMultiFit
    //   GSLSimAn
    //   Genetic
-   ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+   ROOT::Math::Minimizer* min = 
+      ROOT::Math::Factory::CreateMinimizer(minName, algoName);
 
    // set tolerance , etc...
    min->SetMaxFunctionCalls(1000000); // for Minuit/Minuit2 
@@ -70,9 +77,11 @@ int NumericalMinimization(const char * minName = "Minuit2",const char *algoName
 
    // expected minimum is 0
    if ( min->MinValue()  < 1.E-4  && f(xs) < 1.E-4) 
-      std::cout << "Minimizer " << minName << " - " << algoName << "   converged to the right minimum" << std::endl;
+      std::cout << "Minimizer " << minName << " - " << algoName 
+                << "   converged to the right minimum" << std::endl;
    else {
-      std::cout << "Minimizer " << minName << " - " << algoName << "   failed to converge !!!" << std::endl;
+      std::cout << "Minimizer " << minName << " - " << algoName 
+                << "   failed to converge !!!" << std::endl;
       Error("NumericalMinimization","fail to converge");
    }
  

tutorials/fit/combinedFit.C

+//+ Combined (simultaneous) fit of two histogram with separate functions 
+//  and some common parameters
+//
+// See http://root.cern.ch/phpBB3//viewtopic.php?f=3&t=11740#p50908
+// for a modified version working with Fumili or GSLMultiFit 
+//
+// N.B. this macro must be compiled with ACliC 
+//
+//Author: L. Moneta - Dec 2010
+
+#include "Fit/Fitter.h"
+#include "Fit/BinData.h"
+#include "Fit/Chi2FCN.h"
+#include "TH1.h"
+#include "TList.h"
+#include "Math/WrappedMultiTF1.h"
+#include "HFitInterface.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+
+
+// definition of shared parameter
+// background function 
+int iparB[2] = { 0,      // exp amplitude in B histo
+                 2    // exp common parameter 
+};
+
+// signal + background function 
+int iparSB[5] = { 1, // exp amplitude in S+B histo
+                  2, // exp common parameter
+                  3, // gaussian amplitude
+                  4, // gaussian mean
+                  5  // gaussian sigma
+};
+
+struct GlobalChi2 { 
+   GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,  
+                ROOT::Math::IMultiGenFunction & f2) : 
+      fChi2_1(&f1), fChi2_2(&f2) {}
+
+   // parameter vector is first background (in common 1 and 2) 
+   // and then is signal (only in 2)
+   double operator() (const double *par) const {
+      double p1[2];
+      for (int i = 0; i < 2; ++i) p1[i] = par[iparB[i] ];
+
+      double p2[5]; 
+      for (int i = 0; i < 5; ++i) p2[i] = par[iparSB[i] ];
+
+      return (*fChi2_1)(p1) + (*fChi2_2)(p2);
+   } 
+
+   const  ROOT::Math::IMultiGenFunction * fChi2_1;
+   const  ROOT::Math::IMultiGenFunction * fChi2_2;
+};
+
+void combinedFit() { 
+
+  TH1D * hB = new TH1D("hB","histo B",100,0,100);
+  TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100);
+
+  TF1 * fB = new TF1("fB","expo",0,100);
+  fB->SetParameters(1,-0.05);
+  hB->FillRandom("fB");
+
+  TF1 * fS = new TF1("fS","gaus",0,100);
+  fS->SetParameters(1,30,5);
+
+  hSB->FillRandom("fB",2000);
+  hSB->FillRandom("fS",1000);
+
+  // perform now global fit 
+
+  TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100);
+
+  ROOT::Math::WrappedMultiTF1 wfB(*fB,1);
+  ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1);
+
+  ROOT::Fit::DataOptions opt; 
+  ROOT::Fit::DataRange rangeB; 
+  // set the data range
+  rangeB.SetRange(10,90);
+  ROOT::Fit::BinData dataB(opt,rangeB); 
+  ROOT::Fit::FillData(dataB, hB);
+
+  ROOT::Fit::DataRange rangeSB; 
+  rangeSB.SetRange(10,50);
+  ROOT::Fit::BinData dataSB(opt,rangeSB); 
+  ROOT::Fit::FillData(dataSB, hSB);
+
+  ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
+  ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);
+
+  GlobalChi2 globalChi2(chi2_B, chi2_SB);
+
+  ROOT::Fit::Fitter fitter;
+
+  const int Npar = 6; 
+  double par0[Npar] = { 5,5,-0.1,100, 30,10};
+
+  // create before the parameter settings in order to fix or set range on them
+  fitter.Config().SetParamsSettings(6,par0);
+  // fix 5-th parameter  
+  fitter.Config().ParSettings(4).Fix();
+  // set limits on the third and 4-th parameter
+  fitter.Config().ParSettings(2).SetLimits(-10,-1.E-4);
+  fitter.Config().ParSettings(3).SetLimits(0,10000);
+  fitter.Config().ParSettings(3).SetStepSize(5);
+
+  fitter.Config().MinimizerOptions().SetPrintLevel(0);
+  fitter.Config().SetMinimizer("Minuit2","Migrad"); 
+
+  // fit FCN function directly 
+  // (specify optionally data size and flag to indicate that is a chi2 fit)
+  fitter.FitFCN(6,globalChi2,0,dataB.Size()+dataSB.Size(),true);
+  ROOT::Fit::FitResult result = fitter.Result();
+  result.Print(std::cout);
+
+  TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two histograms",
+                             10,10,700,700);
+  c1->Divide(1,2);
+  c1->cd(1);
+  gStyle->SetOptFit(1111);
+
+  fB->SetFitResult( result, iparB);
+  fB->SetRange(rangeB().first, rangeB().second);   
+  fB->SetLineColor(kBlue);
+  hB->GetListOfFunctions()->Add(fB);
+  hB->Draw(); 
+
+  c1->cd(2);
+  fSB->SetFitResult( result, iparSB);
+  fSB->SetRange(rangeSB().first, rangeSB().second);   
+  fSB->SetLineColor(kRed);
+  hSB->GetListOfFunctions()->Add(fSB);
+  hSB->Draw(); 
+    
+
+}