double frac1 =0.2;
double frac2 =0.15;
double frac3 =0.65;

int nbin =40; 
double lo=0, hi=10;

double m1 = 2, m2=3.1, m3 =6;
double s1 = 0.5, s2 = 0.7, s3 =2.7;

void MC(TRandom3* rand, TH1F* h_f1, TH1F* h_f2, TH1F* h_f3){
	int nTotal=100000;
	int nMC =  1000000;

	TH1D* h1 = new TH1D("h1","1st MC", nbin ,lo, hi);
	TH1D* h2 = new TH1D("h2","2nd MC", nbin ,lo, hi);
	TH1D* h3 = new TH1D("h3","3rd MC", nbin ,lo, hi);
	TH1D* hData = new TH1D("htotal","data", nbin ,lo, hi);

	for (int j=0; j<nMC; j++){
		double x=-1;
		do{
			x = rand->Gaus(m1,s1);
		}while(lo>x||x>hi );
		h1->Fill(x);
		do{
			x = rand->Gaus(m2,s2);
		}while(lo>x||x>hi );
		h2->Fill(x);
		do{
			x = rand->Gaus(m3,s3);
		}while(lo>x||x>hi );
		h3->Fill(x);
	}

	for (int j=0; j<nTotal; j++){
		double tmp=rand->Rndm();
		double x =-1;
		if (tmp<frac1){
			do{
				x = rand->Gaus(m1,s1);
			}while(lo>x||x>hi );
			hData->Fill(x);
		}else if (tmp<(frac1+frac2)){
			do{
				x = rand->Gaus(m2,s2);
			}while(lo>x||x>hi );
			hData->Fill(x);
		}else{
			do{
				x = rand->Gaus(m3,s3);
			}while(lo>x||x>hi );
			hData->Fill(x);
		}

	}

	TObjArray *mc = new TObjArray(3);
	mc->Add(h1);
	mc->Add(h2);
	mc->Add(h3);

	TFractionFitter* fit = new TFractionFitter(hData, mc);
	Int_t status = fit->Fit();
	if (status == 0) {                       
		double f1 , f2 , f3;
		double e1, e2, e3;

		fit->GetResult(0,f1,e1);
		fit->GetResult(1,f2,e2);
		fit->GetResult(2,f3,e3);

		h_f1->Fill((frac1-f1)/e1);
		h_f2->Fill((frac2-f2)/e2);
		h_f3->Fill((frac3-f3)/e3);

	}
	delete h1;
	delete h2;
	delete h3;
	delete hData;
	delete mc;
	delete fit;
}

void Draw(){
	int nTotal=10000;
	TRandom3 * rand =new TRandom3();

	TH1D* hData = new TH1D("htotal","data", nbin ,lo, hi);
	TH1D* h1T = new TH1D("h1T","1st MC", nbin ,lo, hi);
	TH1D* h2T = new TH1D("h2T","2nd MC", nbin ,lo, hi);
	TH1D* h3T = new TH1D("h3T","3rd MC", nbin ,lo, hi);
	for (int j=0; j<nTotal; j++){
		double tmp=rand->Rndm();
		double x =-1;
		if (tmp<frac1){
			do{
				x = rand->Gaus(m1,s1);
			}while(lo>x||x>hi );
			hData->Fill(x);
			h1T->Fill(x);
		}else if (tmp<frac1+frac2){
			do{
				x = rand->Gaus(m2,s2);
			}while(lo>x||x>hi );
			hData->Fill(x);
			h2T->Fill(x);
		}else{
			do{
				x = rand->Gaus(m3,s3);
			}while(lo>x||x>hi );
			x = rand->Gaus(m3,s3);
			hData->Fill(x);
			h3T->Fill(x);
		}

	}
	TCanvas* c1 = new TCanvas();
	hData->SetMarkerStyle(21);
	hData->SetMarkerColor(kBlack);
	hData->SetLineColor(kBlack);
	hData->Draw("EP");
	h1T->SetLineColor(kRed);
	h1T->Draw("same");
	h2T->SetLineColor(kOrange+1);
	h2T->Draw("same");
	h3T->SetLineColor(kBlue);
	h3T->Draw("same");
	c1->SaveAs("data.png");


}

void toyMC(){
	int nTrial =10000;

	TH1F* h_f1 = new TH1F("h_f1","frac1;(true-fit)/error",40,-1.5,1.5);
	TH1F* h_f2 = new TH1F("h_f2","frac2;(true-fit)/error",40,-1.5,1.5);
	TH1F* h_f3 = new TH1F("h_f3","frac3;(true-fit)/error",40,-1.5,1.5);

	TRandom3 * rand =new TRandom3();

	for (int i=0; i<nTrial; i++){
		MC(rand, h_f1, h_f2, h_f3);
	}

	gStyle->SetOptFit(1);
	gStyle->SetOptStat(0);
	TCanvas* c1 = new TCanvas();
	h_f1->Draw();
	h_f1->Fit("gaus","L");
	c1->SaveAs("frac1_error.png");
	TCanvas* c2 = new TCanvas();
	h_f2->Draw();
	h_f2->Fit("gaus","L");
	c2->SaveAs("frac2_error.png");
	TCanvas* c3 = new TCanvas();
	h_f3->Draw();
	h_f3->Fit("gaus","L");
	c3->SaveAs("frac3_error.png");


	Draw();
}