ROOT logo

From $ROOTSYS/tutorials/fit/TwoHistoFit2D.C

//+ Example to fit two histograms at the same time 
//Author: Rene Brun

#include "TH2D.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "TVirtualFitter.h"
#include "TList.h"

#include <vector>
#include <map>
#include <iostream>
  
double gauss2D(double *x, double *par) {
   double z1 = double((x[0]-par[1])/par[2]);
   double z2 = double((x[1]-par[3])/par[4]);
   return par[0]*exp(-0.5*(z1*z1+z2*z2));
}   
double my2Dfunc(double *x, double *par) {
   double *p1 = &par[0];
   double *p2 = &par[5];
   return gauss2D(x,p1) + gauss2D(x,p2);
}



// data need to be globals to be visible by fcn 

std::vector<std::pair<double, double> > coords;        
std::vector<double > values;        
std::vector<double > errors;        

void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */  )
{
  int n = coords.size();
  double chi2 = 0; 
  double tmp,x[2]; 
  for (int i = 0; i <n; ++i ) { 
    x[0] = coords[i].first;
    x[1] = coords[i].second;
    tmp = ( values[i] - my2Dfunc(x,p))/errors[i];
    chi2 += tmp*tmp; 
  }
  fval = chi2; 
}
TRandom3 rndm; 
void FillHisto(TH2D * h, int n, double * p) { 


  const double mx1 = p[1]; 
  const double my1 = p[3]; 
  const double sx1 = p[2]; 
  const double sy1 = p[4]; 
  const double mx2 = p[6]; 
  const double my2 = p[8]; 
  const double sx2 = p[7]; 
  const double sy2 = p[9]; 
  //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2); 
  const double w1 = 0.5; 

  double x, y; 
  for (int i = 0; i < n; ++i) {
    // generate randoms with larger gaussians
    rndm.Rannor(x,y);

    double r = rndm.Rndm(1);
    if (r < w1) { 
      x = x*sx1 + mx1; 
      y = y*sy1 + my1; 
    }
    else { 
      x = x*sx2 + mx2; 
      y = y*sy2 + my2; 
    }      
    h->Fill(x,y);
     
  }
}
      
  


int TwoHistoFit2D(bool global = true) { 

  // create two histograms 

  int nbx1 = 50;
  int nby1 = 50;
  int nbx2 = 50;
  int nby2 = 50;
  double xlow1 = 0.; 
  double ylow1 = 0.; 
  double xup1 = 10.; 
  double yup1 = 10.; 
  double xlow2 = 5.; 
  double ylow2 = 5.; 
  double xup2 = 20.; 
  double yup2 = 20.; 

  TH2D * h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
  TH2D * h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2);

  double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. };
  // create fit function
  TF2 * func = new TF2("func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);
  func->SetParameters(iniParams);

  // fill Histos
  int n1 = 1000000;
  int n2 = 1000000; 
  //  h1->FillRandom("func", n1);
  //h2->FillRandom("func",n2);
  FillHisto(h1,n1,iniParams);
  FillHisto(h2,n2,iniParams);

  // scale histograms to same heights (for fitting)
  double dx1 = (xup1-xlow1)/double(nbx1); 
  double dy1 = (yup1-ylow1)/double(nby1);
  double dx2 = (xup2-xlow2)/double(nbx2);
  double dy2 = (yup2-ylow2)/double(nby2);
//   h1->Sumw2();
//   h1->Scale( 1.0 / ( n1 * dx1 * dy1 ) );
  // scale histo 2 to scale of 1 
  h2->Sumw2();
  h2->Scale(  ( double(n1) * dx1 * dy1 )  / ( double(n2) * dx2 * dy2 ) );


  if (global) { 
    // fill data structure for fit (coordinates + values + errors) 
    std::cout << "Do global fit" << std::endl;
    // fit now all the function together
    
    // fill data structure for fit (coordinates + values + errors) 
    TAxis *xaxis1  = h1->GetXaxis();
    TAxis *yaxis1  = h1->GetYaxis();
    TAxis *xaxis2  = h2->GetXaxis();
    TAxis *yaxis2  = h2->GetYaxis();

    int nbinX1 = h1->GetNbinsX(); 
    int nbinY1 = h1->GetNbinsY(); 
    int nbinX2 = h2->GetNbinsX(); 
    int nbinY2 = h2->GetNbinsY(); 

  /// reset data structure
  coords = std::vector<std::pair<double,double> >();
  values = std::vector<double>();
  errors = std::vector<double>();
  

  for (int ix = 1; ix <= nbinX1; ++ix) { 
    for (int iy = 1; iy <= nbinY1; ++iy) { 
      if ( h1->GetBinContent(ix,iy) > 0 ) { 
         coords.push_back( std::make_pair(xaxis1->GetBinCenter(ix), yaxis1->GetBinCenter(iy) ) );
         values.push_back( h1->GetBinContent(ix,iy) );
         errors.push_back( h1->GetBinError(ix,iy) );
      }
    }
  }
  for (int ix = 1; ix <= nbinX2; ++ix) { 
    for (int iy = 1; iy <= nbinY2; ++iy) { 
      if ( h2->GetBinContent(ix,iy) > 0 ) { 
         coords.push_back( std::make_pair(xaxis2->GetBinCenter(ix), yaxis2->GetBinCenter(iy) ) );
         values.push_back( h2->GetBinContent(ix,iy) );
         errors.push_back( h2->GetBinError(ix,iy) );
      }
    }
  }

  TVirtualFitter::SetDefaultFitter("Minuit");
  TVirtualFitter * minuit = TVirtualFitter::Fitter(0,10);
  for (int i = 0; i < 10; ++i) {  
    minuit->SetParameter(i, func->GetParName(i), func->GetParameter(i), 0.01, 0,0);
  }
  minuit->SetFCN(myFcn);

  double arglist[100];
  arglist[0] = 0;
  // set print level
  minuit->ExecuteCommand("SET PRINT",arglist,2);

 // minimize
  arglist[0] = 5000; // number of function calls
  arglist[1] = 0.01; // tolerance
  minuit->ExecuteCommand("MIGRAD",arglist,2);

  //get result
  double minParams[10];
  double parErrors[10];
  for (int i = 0; i < 10; ++i) {  
    minParams[i] = minuit->GetParameter(i);
    parErrors[i] = minuit->GetParError(i);
  }
  double chi2, edm, errdef; 
  int nvpar, nparx;
  minuit->GetStats(chi2,edm,errdef,nvpar,nparx);

  func->SetParameters(minParams);
  func->SetParErrors(parErrors);
  func->SetChisquare(chi2);
  int ndf = coords.size()-nvpar;
  func->SetNDF(ndf);

  std::cout << "Chi2 Fit = " << chi2 << " ndf = " << ndf << "  " << func->GetNDF() << std::endl;

  // add to list of functions
  h1->GetListOfFunctions()->Add(func);
  h2->GetListOfFunctions()->Add(func);
  }
  else {     
    // fit independently
    h1->Fit(func);
    h2->Fit(func);
  }

     

  // Create a new canvas.
  TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800);
  c1->Divide(2,2);
  gStyle->SetOptFit();
  gStyle->SetStatY(0.6);

  c1->cd(1);
  h1->Draw();
  func->SetRange(xlow1,ylow1,xup1,yup1);
  func->DrawCopy("cont1 same");
  c1->cd(2);
  h1->Draw("lego");
  func->DrawCopy("surf1 same");
  c1->cd(3);
  func->SetRange(xlow2,ylow2,xup2,yup2);
  h2->Draw();
  func->DrawCopy("cont1 same");
  c1->cd(4);
  h2->Draw("lego");
  gPad->SetLogz();
  func->Draw("surf1 same");

  return 0; 
}
 TwoHistoFit2D.C:1
 TwoHistoFit2D.C:2
 TwoHistoFit2D.C:3
 TwoHistoFit2D.C:4
 TwoHistoFit2D.C:5
 TwoHistoFit2D.C:6
 TwoHistoFit2D.C:7
 TwoHistoFit2D.C:8
 TwoHistoFit2D.C:9
 TwoHistoFit2D.C:10
 TwoHistoFit2D.C:11
 TwoHistoFit2D.C:12
 TwoHistoFit2D.C:13
 TwoHistoFit2D.C:14
 TwoHistoFit2D.C:15
 TwoHistoFit2D.C:16
 TwoHistoFit2D.C:17
 TwoHistoFit2D.C:18
 TwoHistoFit2D.C:19
 TwoHistoFit2D.C:20
 TwoHistoFit2D.C:21
 TwoHistoFit2D.C:22
 TwoHistoFit2D.C:23
 TwoHistoFit2D.C:24
 TwoHistoFit2D.C:25
 TwoHistoFit2D.C:26
 TwoHistoFit2D.C:27
 TwoHistoFit2D.C:28
 TwoHistoFit2D.C:29
 TwoHistoFit2D.C:30
 TwoHistoFit2D.C:31
 TwoHistoFit2D.C:32
 TwoHistoFit2D.C:33
 TwoHistoFit2D.C:34
 TwoHistoFit2D.C:35
 TwoHistoFit2D.C:36
 TwoHistoFit2D.C:37
 TwoHistoFit2D.C:38
 TwoHistoFit2D.C:39
 TwoHistoFit2D.C:40
 TwoHistoFit2D.C:41
 TwoHistoFit2D.C:42
 TwoHistoFit2D.C:43
 TwoHistoFit2D.C:44
 TwoHistoFit2D.C:45
 TwoHistoFit2D.C:46
 TwoHistoFit2D.C:47
 TwoHistoFit2D.C:48
 TwoHistoFit2D.C:49
 TwoHistoFit2D.C:50
 TwoHistoFit2D.C:51
 TwoHistoFit2D.C:52
 TwoHistoFit2D.C:53
 TwoHistoFit2D.C:54
 TwoHistoFit2D.C:55
 TwoHistoFit2D.C:56
 TwoHistoFit2D.C:57
 TwoHistoFit2D.C:58
 TwoHistoFit2D.C:59
 TwoHistoFit2D.C:60
 TwoHistoFit2D.C:61
 TwoHistoFit2D.C:62
 TwoHistoFit2D.C:63
 TwoHistoFit2D.C:64
 TwoHistoFit2D.C:65
 TwoHistoFit2D.C:66
 TwoHistoFit2D.C:67
 TwoHistoFit2D.C:68
 TwoHistoFit2D.C:69
 TwoHistoFit2D.C:70
 TwoHistoFit2D.C:71
 TwoHistoFit2D.C:72
 TwoHistoFit2D.C:73
 TwoHistoFit2D.C:74
 TwoHistoFit2D.C:75
 TwoHistoFit2D.C:76
 TwoHistoFit2D.C:77
 TwoHistoFit2D.C:78
 TwoHistoFit2D.C:79
 TwoHistoFit2D.C:80
 TwoHistoFit2D.C:81
 TwoHistoFit2D.C:82
 TwoHistoFit2D.C:83
 TwoHistoFit2D.C:84
 TwoHistoFit2D.C:85
 TwoHistoFit2D.C:86
 TwoHistoFit2D.C:87
 TwoHistoFit2D.C:88
 TwoHistoFit2D.C:89
 TwoHistoFit2D.C:90
 TwoHistoFit2D.C:91
 TwoHistoFit2D.C:92
 TwoHistoFit2D.C:93
 TwoHistoFit2D.C:94
 TwoHistoFit2D.C:95
 TwoHistoFit2D.C:96
 TwoHistoFit2D.C:97
 TwoHistoFit2D.C:98
 TwoHistoFit2D.C:99
 TwoHistoFit2D.C:100
 TwoHistoFit2D.C:101
 TwoHistoFit2D.C:102
 TwoHistoFit2D.C:103
 TwoHistoFit2D.C:104
 TwoHistoFit2D.C:105
 TwoHistoFit2D.C:106
 TwoHistoFit2D.C:107
 TwoHistoFit2D.C:108
 TwoHistoFit2D.C:109
 TwoHistoFit2D.C:110
 TwoHistoFit2D.C:111
 TwoHistoFit2D.C:112
 TwoHistoFit2D.C:113
 TwoHistoFit2D.C:114
 TwoHistoFit2D.C:115
 TwoHistoFit2D.C:116
 TwoHistoFit2D.C:117
 TwoHistoFit2D.C:118
 TwoHistoFit2D.C:119
 TwoHistoFit2D.C:120
 TwoHistoFit2D.C:121
 TwoHistoFit2D.C:122
 TwoHistoFit2D.C:123
 TwoHistoFit2D.C:124
 TwoHistoFit2D.C:125
 TwoHistoFit2D.C:126
 TwoHistoFit2D.C:127
 TwoHistoFit2D.C:128
 TwoHistoFit2D.C:129
 TwoHistoFit2D.C:130
 TwoHistoFit2D.C:131
 TwoHistoFit2D.C:132
 TwoHistoFit2D.C:133
 TwoHistoFit2D.C:134
 TwoHistoFit2D.C:135
 TwoHistoFit2D.C:136
 TwoHistoFit2D.C:137
 TwoHistoFit2D.C:138
 TwoHistoFit2D.C:139
 TwoHistoFit2D.C:140
 TwoHistoFit2D.C:141
 TwoHistoFit2D.C:142
 TwoHistoFit2D.C:143
 TwoHistoFit2D.C:144
 TwoHistoFit2D.C:145
 TwoHistoFit2D.C:146
 TwoHistoFit2D.C:147
 TwoHistoFit2D.C:148
 TwoHistoFit2D.C:149
 TwoHistoFit2D.C:150
 TwoHistoFit2D.C:151
 TwoHistoFit2D.C:152
 TwoHistoFit2D.C:153
 TwoHistoFit2D.C:154
 TwoHistoFit2D.C:155
 TwoHistoFit2D.C:156
 TwoHistoFit2D.C:157
 TwoHistoFit2D.C:158
 TwoHistoFit2D.C:159
 TwoHistoFit2D.C:160
 TwoHistoFit2D.C:161
 TwoHistoFit2D.C:162
 TwoHistoFit2D.C:163
 TwoHistoFit2D.C:164
 TwoHistoFit2D.C:165
 TwoHistoFit2D.C:166
 TwoHistoFit2D.C:167
 TwoHistoFit2D.C:168
 TwoHistoFit2D.C:169
 TwoHistoFit2D.C:170
 TwoHistoFit2D.C:171
 TwoHistoFit2D.C:172
 TwoHistoFit2D.C:173
 TwoHistoFit2D.C:174
 TwoHistoFit2D.C:175
 TwoHistoFit2D.C:176
 TwoHistoFit2D.C:177
 TwoHistoFit2D.C:178
 TwoHistoFit2D.C:179
 TwoHistoFit2D.C:180
 TwoHistoFit2D.C:181
 TwoHistoFit2D.C:182
 TwoHistoFit2D.C:183
 TwoHistoFit2D.C:184
 TwoHistoFit2D.C:185
 TwoHistoFit2D.C:186
 TwoHistoFit2D.C:187
 TwoHistoFit2D.C:188
 TwoHistoFit2D.C:189
 TwoHistoFit2D.C:190
 TwoHistoFit2D.C:191
 TwoHistoFit2D.C:192
 TwoHistoFit2D.C:193
 TwoHistoFit2D.C:194
 TwoHistoFit2D.C:195
 TwoHistoFit2D.C:196
 TwoHistoFit2D.C:197
 TwoHistoFit2D.C:198
 TwoHistoFit2D.C:199
 TwoHistoFit2D.C:200
 TwoHistoFit2D.C:201
 TwoHistoFit2D.C:202
 TwoHistoFit2D.C:203
 TwoHistoFit2D.C:204
 TwoHistoFit2D.C:205
 TwoHistoFit2D.C:206
 TwoHistoFit2D.C:207
 TwoHistoFit2D.C:208
 TwoHistoFit2D.C:209
 TwoHistoFit2D.C:210
 TwoHistoFit2D.C:211
 TwoHistoFit2D.C:212
 TwoHistoFit2D.C:213
 TwoHistoFit2D.C:214
 TwoHistoFit2D.C:215
 TwoHistoFit2D.C:216
 TwoHistoFit2D.C:217
 TwoHistoFit2D.C:218
 TwoHistoFit2D.C:219
 TwoHistoFit2D.C:220
 TwoHistoFit2D.C:221
 TwoHistoFit2D.C:222
 TwoHistoFit2D.C:223
 TwoHistoFit2D.C:224
 TwoHistoFit2D.C:225
 TwoHistoFit2D.C:226
 TwoHistoFit2D.C:227
 TwoHistoFit2D.C:228
 TwoHistoFit2D.C:229
 TwoHistoFit2D.C:230
 TwoHistoFit2D.C:231
 TwoHistoFit2D.C:232
 TwoHistoFit2D.C:233
 TwoHistoFit2D.C:234
 TwoHistoFit2D.C:235
 TwoHistoFit2D.C:236
 TwoHistoFit2D.C:237
 TwoHistoFit2D.C:238
 TwoHistoFit2D.C:239
 TwoHistoFit2D.C:240
 TwoHistoFit2D.C:241
 TwoHistoFit2D.C:242
 TwoHistoFit2D.C:243