// pointers to the data
TH1F *data;                              //data histogram
TH1F *mc0;                               // first MC histogram
TH1F *mc1;                               // second MC histogram
TH1F *mc2;                               // third MC histogram

// Global canvas and histos for display
TCanvas *c;
TH1F *hff0, *hc20, *hlh0;
TH1F *hff1, *hc21, *hlh1;
TH1F *hff2, *hc22, *hlh2;

//------------------------------------------------------------------------------

void chi2(Int_t &npar, Double_t *gin, Double_t &fout, Double_t *par, Int_t iflag) {	
	// simple chi2
	// par[] is the array of the parameters = proportions
	// fout is the output value

	Double_t ni, s0i, s1i, s2i, fi;
	Double_t N  = data->GetEntries();
	Double_t N0 = mc0->GetEntries();
	Double_t N1 = mc1->GetEntries();
	Double_t N2 = mc2->GetEntries();
	
	fout = 0.;
	for( Short_t iBin=1; iBin<=data->GetNbinsX(); iBin++ ) {
		ni  = data->GetBinContent( iBin);
		s0i = mc0->GetBinContent( iBin);
		s1i = mc1->GetBinContent( iBin);
		s2i = mc2->GetBinContent( iBin);
		fi  = N * ( s0i/N0*par[0] + s1i/N1*par[1] + s2i/N2*par[2] );
		if( ni>0. ) fout += (ni-fi)*(ni-fi)/ni;
	}
	
}

//------------------------------------------------------------------------------

void likelyhood(Int_t &npar, Double_t *gin, Double_t &fout, Double_t *par, Int_t iflag) {	
	// simple ln(likelyhood)
	// par[] is the array of the parameters = proportions
	// fout is the output value
	
	Double_t ni, s0i, s1i, s2i, fi;
	Double_t N  = data->GetEntries();
	Double_t N0 = mc0->GetEntries();
	Double_t N1 = mc1->GetEntries();
	Double_t N2 = mc2->GetEntries();
	
	fout = 0.;
	for( Short_t iBin=1; iBin<=data->GetNbinsX(); iBin++ ) {
		ni  = data->GetBinContent( iBin);
		s0i = mc0->GetBinContent( iBin);
		s1i = mc1->GetBinContent( iBin);
		s2i = mc2->GetBinContent( iBin);
		fi  = N * ( s0i/N0*par[0] + s1i/N1*par[1] + s2i/N2*par[2] );
		fout -= ni*TMath::Log(fi)-fi;
	}
	
}

//------------------------------------------------------------------------------
void fractionFitterComparison(){
	
	// Comparing TFractionFitter class with simpler chi2 and likelyhood
	// 1 Dimension only, x is an angle fron 0 to pi
	//
	// JB May 2010
	
	Int_t nFits = 1000; // number of fits to perform
	
	// parameters and functions to generate the data
	Int_t Ndata = 1000;
	Int_t N0 = 10000;
	Int_t N1 = 10000;
	Int_t N2 = 10000;
	
	Int_t nBins = 40;
	
	Double_t trueP0 = .1;
	Double_t trueP1 = .3;
	Double_t trueP2 = 1.-trueP0-trueP1;
	
	// Histo of solution for p0
	hff0 = new TH1F("hff0", "Solution from FractionFitter", 50, 0, 1.);
	hff0->SetXTitle("p0");
	hc20 = new TH1F("hc20", "Solution from simple #chi 2", 50, 0, 1.);
	hc20->SetXTitle("p0");
	hlh0 = new TH1F("hlh0", "Solution from simple likelyhood", 50, 0, 1.);
	hlh0->SetXTitle("p0");

	hff1 = new TH1F("hff1", "Solution from FractionFitter", 50, 0, 1.);
	hff1->SetXTitle("p1");
	hc21 = new TH1F("hc21", "Solution from simple #chi 2", 50, 0, 1.);
	hc21->SetXTitle("p1");
	hlh1 = new TH1F("hlh1", "Solution from simple likelyhood", 50, 0, 1.);
	hlh1->SetXTitle("p1");

	hff2 = new TH1F("hff2", "Solution from FractionFitter", 50, 0, 1.);
	hff2->SetXTitle("p2");
	hc22 = new TH1F("hc22", "Solution from simple #chi 2", 50, 0, 1.);
	hc22->SetXTitle("p2");
	hlh2 = new TH1F("hlh2", "Solution from simple likelyhood", 50, 0, 1.);
	hlh2->SetXTitle("p2");
	
	// contribution 0
	TF1 *f0 = new TF1("f0", "[0]*(1-cos(x))/TMath::Pi()", 0., TMath::Pi());
	f0->SetParameter(0,1.);
	f0->SetLineColor(2);
	Double_t int0 = f0->Integral( 0., TMath::Pi());

	// contribution 1
	TF1 *f1 = new TF1("f1", "[0]*(1-cos(x)*cos(x))*2./TMath::Pi()", 0., TMath::Pi());
	f1->SetParameter(0,1.);
	f1->SetLineColor(3);
	Double_t int1 = f1->Integral( 0., TMath::Pi());

	// contribution 2
	TF1 *f2 = new TF1("f2", "[0]*(1+cos(x))/TMath::Pi()", 0., TMath::Pi());
	f2->SetParameter(0,1.);
	f2->SetLineColor(4);
	Double_t int2 = f2->Integral( 0., TMath::Pi());

	// ***************** FIT LOOP
	Double_t p, x;	
	TFractionFitter* fit;
	TFitter *minimizer;
	for( Int_t iFit=0; iFit<nFits; iFit++) {
		
	
	// generate data
	data = new TH1F("data", "Data angle distribution", nBins, 0, TMath::Pi());
	data->SetXTitle("x");
	for( Int_t i=0; i<Ndata; i++) {
		p = gRandom->Uniform();
		if( p<trueP0 ) { 
			x = f0->GetRandom();
		}
		else if( p<trueP0+trueP1 ) { 
			x = f1->GetRandom();
		}
		else { 
			x = f2->GetRandom();
		}
		data->Fill( x);
	}
	
	
	// generate MC samples
	mc0 = new TH1F("mc0", "MC sample 0 angle distribution", nBins, 0, TMath::Pi());
	mc0->SetXTitle("x");
	for( Int_t i=0; i<N0; i++) {
		mc0->Fill( f0->GetRandom() ); 
	}

	mc1 = new TH1F("mc1", "MC sample 1 angle distribution", nBins, 0, TMath::Pi());
	mc1->SetXTitle("x");
	for( Int_t i=0; i<N1; i++) {
		mc1->Fill( f1->GetRandom() ); 
	}

	mc2 = new TH1F("mc2", "MC sample 2 angle distribution", nBins, 0, TMath::Pi());
	mc2->SetXTitle("x");
	for( Int_t i=0; i<N2; i++) {
		mc2->Fill( f2->GetRandom() ); 
	}

	cout << endl << "Data & MC generation DONE" << endl << endl;
	
	// FractionFitter
	cout << endl << "Using FractionFitter " << endl;
	TObjArray *mc = new TObjArray(3);        // MC histograms are put in this array
	mc->Add(mc0);
	mc->Add(mc1);
	mc->Add(mc2);
	fit = new TFractionFitter(data, mc); // initialise
	fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
	fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
	fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
	//fit->SetRangeX(1,15);                    // use only the first 15 bins in the fit
	Int_t status = fit->Fit();               // perform the fit
	cout << "FractionFitter fit status: " << status << endl;
	Double_t ffp0, ffp1, ffp2, fferrP0, fferrP1, fferrP2;
	if (status == 0) {                       // check on fit status
		TH1F* result = (TH1F*) fit->GetPlot();
		fit->GetResult( 0, ffp0, fferrP0);
		fit->GetResult( 1, ffp1, fferrP1);
		fit->GetResult( 2, ffp2, fferrP2);
		printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 0, trueP0, ffp0, fferrP0);
		printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 1, trueP1, ffp1, fferrP1);
		printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 2, trueP2, ffp2, fferrP2);
		hff0->Fill( ffp0);
		hff1->Fill( ffp1);
		hff2->Fill( ffp2);
	}
	delete mc;
	delete fit;
		
	// chi2
	cout << endl << "Using Simple chi2 " << endl;
	minimizer = new TFitter(3);
	minimizer->SetFCN( chi2); // function to minimize
	minimizer->SetParameter(0,"p0",.5,.1,0,0);
	minimizer->SetParameter(1,"p1",.5,.1,0,0);
	minimizer->SetParameter(2,"p2",.5,.1,0,0);
	minimizer->ExecuteCommand("MIGRAD",0,0);
	Double_t  c2p0, c2p1, c2p2, c2errP0, c2errP1, c2errP2;
	c2p0 = minimizer->GetParameter(0);
	c2errP0 = minimizer->GetParError(0);
	c2p1 = minimizer->GetParameter(1);
	c2errP1 = minimizer->GetParError(1);
	c2p2 = minimizer->GetParameter(2);
	c2errP2 = minimizer->GetParError(2);
	cout << "Simple chi2 fit status: " << endl;
	printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 0, trueP0, c2p0, c2errP0);
	printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 1, trueP1, c2p1, c2errP1);
	printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 2, trueP2, c2p2, c2errP2);
	hc20->Fill( c2p0);
	hc21->Fill( c2p1);
	hc22->Fill( c2p2);
	delete minimizer;
		
		// likelyhood
  	cout << endl << "Using Simple likelyhood " << endl;
		minimizer = new TFitter(3);
		minimizer->SetFCN( likelyhood); // function to minimize
		minimizer->SetParameter(0,"p0",.5,.1,0,0);
		minimizer->SetParameter(1,"p1",.5,.1,0,0);
		minimizer->SetParameter(2,"p2",.5,.1,0,0);
    minimizer->ExecuteCommand("SIMPLEX",0,0);
		minimizer->ExecuteCommand("MIGRAD",0,0);
		Double_t  lhp0, lhp1, lhp2, lherrP0, lherrP1, lherrP2;
		lhp0 = minimizer->GetParameter(0);
		lherrP0 = minimizer->GetParError(0);
		lhp1 = minimizer->GetParameter(1);
		lherrP1 = minimizer->GetParError(1);
		lhp2 = minimizer->GetParameter(2);
		lherrP2 = minimizer->GetParError(2);
	cout << "Simple likelyhood fit status: " << endl;
	printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 0, trueP0, lhp0, lherrP0);
	printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 1, trueP1, lhp1, lherrP1);
	printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 2, trueP2, lhp2, lherrP2);
	hlh0->Fill( lhp0);
	hlh1->Fill( lhp1);
	hlh2->Fill( lhp2);

		
	// clean a bit
		delete data;
		delete mc0;
		delete mc1;
		delete mc2;
	} // end fits loop
		
		// Display
	c = new TCanvas("c", "FractionFitter comparison", 900, 700);
	c->Divide(3,3);
	
	c->cd(1);
	hff0->Draw();
	c->cd(4);
	hc20->Draw();
	c->cd(7);
	hlh0->Draw();

	c->cd(2);
	hff1->Draw();
	c->cd(5);
	hc21->Draw();
	c->cd(8);
	hlh1->Draw();
	

	c->cd(3);
	hff2->Draw();
	c->cd(6);
	hc22->Draw();
	c->cd(9);
	hlh2->Draw();
	
}