//g++ `root-config --cflags` fit.cpp `root-config --libs` -o fit
using namespace std;
#include <iostream>
#include "TCanvas.h"
#include "TFile.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TFractionFitter.h"
#include "Fit/Fitter.h"
#include "Fit/FitConfig.h"
#include "TF1.h"
#include "TH1.h"
#include "TStyle.h"

int main(void){
  //for generating random numbers
  TRandom3 r;
  //pointers to histos
  TH1F *data;
  TH1F *mc0;
  TH1F *mc1;
  TH1F *mc;
  
  //parameters and functions to generate the data
  Int_t N0=1000000;
  Int_t N1=1000000;

  Int_t nBins=50;
  
  //contribution 0
  TF1 *f0=new TF1("f0","gaus(0)",0,10);
  f0->SetParameter(0,1.);
  f0->SetParameter(1,1.0);
  f0->SetParameter(2,2.0);
  f0->SetLineColor(kRed);
  Double_t int0=f0->Integral(0,10);

  //contribution 1
  TF1 *f1=new TF1("f1","gaus(0)",0,10);
  f1->SetParameter(0,1.);
  f1->SetParameter(1,9.0);
  f1->SetParameter(2,2.0);
  f1->SetLineColor(kGreen);
  Double_t int1=f1->Integral(0,10);
  
  //summing all MC
  mc=new TH1F("mc","mc",nBins,0,10);
  
  //generating mc samples
  mc0=new TH1F("mc0","mc0",nBins,0,10);
  mc0->SetLineColor(kRed);
  for(Int_t i=0;i<N0;i++){
    mc0->Fill(f0->GetRandom());
  }
  mc1=new TH1F("mc1","mc1",nBins,0,10);
  mc1->SetLineColor(kGreen);
  for(Int_t i=0;i<N1;i++){
    mc1->Fill(f1->GetRandom());
  }

  //adding mc samples
  
  //does not work with or without scaling
  //mc0->Scale(0.2);
  //  mc1->Scale(0.8);
  mc->Add(mc0);
  mc->Add(mc1);
  //generate data
  data=new TH1F("data","data",nBins,0,10);
  Double_t p,x;
  TH1F* pull0=new TH1F ("pull0","pull0",nBins,-5,5);
  TH1F* pull1=new TH1F ("pull1","pull1",nBins,-5,5);
  double mc_sig=mc0->Integral()/(mc0->Integral()+mc1->Integral());
  double mc_bgd=mc1->Integral()/(mc0->Integral()+mc1->Integral());
  int number=1000;
  for(int j=0;j<number;j++){
    //int n=r.Poisson(mc->Integral());
    int n=mc->Integral()/100.;
    for(Int_t i=0;i<n;i++){
      data->Fill(mc->GetRandom());
    }
    
    //TFractionFitter
    TObjArray *MC = new TObjArray(2);
    MC->Add(mc0);
    MC->Add(mc1);
    TFractionFitter *fit =new TFractionFitter(data,MC,"Q");
    Int_t status=fit->Fit();
    double error;
    double value;
    if(status==0){
      fit->GetResult(0,value,error);
      pull0->Fill((mc_sig-value)/error);
      fit->GetResult(1,value,error);
      pull1->Fill((mc_bgd-value)/error);
    }
    if(j==number-1){ 
      TCanvas c("c", "FractionFitter example", 700, 700);
      c.Divide(2,2);
      c.cd(1);
      f0->DrawClone();
      f1->DrawClone("same");
      
      c.cd(2);
      mc0->Draw("PE");
      mc1->Draw("PE same");
      mc->Draw("PE same");
      
      c.cd(3);
      pull0->Draw();
      pull0->Fit("gaus");
      gStyle->SetOptFit(1);
      c.cd(4);
      pull1->Draw();
      pull1->Fit("gaus");
      gStyle->SetOptFit(1);
      c.SaveAs("Try.png");
    }
    fit->Delete();
    data->Reset();
  }
}