TNtupleD *nt = 0;

double model(double x, double *par)
{
    double bprime = par[0];
    double fs     = par[1]/nt->GetEntries();
    double mu     = par[2];
    double sigma  = par[3];
    
    double norm  = 1./sqrt(2.*TMath::Pi())/sigma;
    double G     = norm*exp(-0.5 * pow((x-mu)/sigma,2));
    
    return fs * G + (1.-fs) * (1 + bprime*x)/(2. + 2.*bprime);
}

void fcn(int &npar, double *gin, double &f, double *par, int iflag)
{
    f = 0.;
    for (int i=0;i<nt->GetEntries()/100;i++) {
        nt->GetEntry(i);
        double *mass = nt->GetArgs();
        
        double L = model(*mass,par);
        if (L>0.) f -= 2.*log(L);
        else { f = HUGE; return; }
    }
}

void example_06b()
{
    TFile *fin = new TFile("example_data.root");
    nt = (TNtupleD *)fin->Get("nt");
    
    TH1D *hist2 = new TH1D("hist2","Binned data",100,0.,2.);
    nt->Project("hist2","mass","","",nt->GetEntries()/100);
    hist2->SetFillColor(kRed-9);
    hist2->SetStats(false);
    hist2->GetXaxis()->SetTitle("Mass");
    
    TMinuit *gMinuit = new TMinuit(4);
    gMinuit->SetFCN(fcn);
    
    gMinuit->DefineParameter(0, "bprime",-0.3, 1.,  -10.,   10.);
    gMinuit->DefineParameter(1, "area", 2000., 1.,    0.,20000.);
    gMinuit->DefineParameter(2, "mean",  1.00, 1.,   0.5,   1.5);
    gMinuit->DefineParameter(3, "width", 0.05, 1., 0.001,  0.15);

    gMinuit->Command("MIGRAD");
    gMinuit->Command("MIGRAD");
    gMinuit->Command("MINOS");
    
    double par[4],err[4];
    for(int i=0;i<4;i++) gMinuit->GetParameter(i,par[i],err[i]);
    
    TH1F* curve = new TH1F("curve","curve",hist2->GetNbinsX()*5,
                           hist2->GetXaxis()->GetXmin(),hist2->GetXaxis()->GetXmax());
    
    for(int i=1;i<=curve->GetNbinsX();i++) {
        double x = curve->GetBinCenter(i);
        double f = model(x,par);
        double BinWidth = hist2->GetBinWidth(1);
        curve->SetBinContent(i,f*nt->GetEntries()/100*BinWidth);
    }
    curve->SetLineWidth(3);
    
    hist2->Draw();
    curve->Draw("csame");
}