Double_t* dgauss(TH1* h, int low1, int high1, int low2, int high2, Double_t min = -1, Double_t max = -1);

Double_t* dgauss(TH1* h, int low1, int high1, int low2, int high2, Double_t min, Double_t max)
{
     // fits two gaussians to a histogram
     // low1, high1 -> approximate limits of first (lowest peak) gaussian (for parameter estimation)
     // low2, high2 " "
     // returns the array of parameters and their errors (6 paramters followed by 6 errors)

     Double_t* par = new Double_t[6];

     Int_t max_bin = h->GetMaximumBin();
     Int_t min_bin = h->GetMinimumBin();
     
     if(max==-1){
	  max = h->GetBinCenter(max_bin);
     }
     if(min==-1){
	  min = h->GetBinCenter(min_bin);
     }

     TString n1("g1_"); n1+=h->GetName();
     TString n2("g2_"); n2+=h->GetName();
     TString ntot("total_"); ntot+=h->GetName();
     
     TF1* g1 = new TF1(n1.Data(), "gaus", low1, high1);
     TF1* g2 = new TF1(n2.Data(), "gaus", low2, high2);
     TF1* tot = new TF1(ntot.Data(), "gaus(0)+gaus(3)", min, max);

     tot->SetLineColor(1);
     tot->SetLineWidth(3);
     tot->SetParName(0, "Constant 1");
     tot->SetParName(1, "Mean 1");
     tot->SetParName(2, "Sigma 1");
     tot->SetParName(3, "Constant 2");
     tot->SetParName(4, "Mean 2");
     tot->SetParName(5, "Sigma 2");


     g1->SetLineStyle(2);
     g2->SetLineStyle(2);
     g1->SetFillColor(2);
     g2->SetFillColor(4);
     g1->SetFillStyle(3004);
     g2->SetFillStyle(3005);

     // fit for estimation
     gStyle->SetOptStat(0);
     gStyle->SetOptFit(0);

     h->Fit(n1.Data(),"IR0");
     h->Fit(n2.Data(),"IR0+");
     g1->GetParameters(&par[0]);
     g2->GetParameters(&par[3]);
     tot->SetParameters(par);

     g1->SetRange(min,max);
     g2->SetRange(min,max);

     h->Fit(ntot.Data(), "IR+");

     
     tot->GetParameters(par);

     g1->SetParameters(&par[0]);
     g2->SetParameters(&par[3]);
     
     
     h->Draw("e1");
     tot->Draw("csame");

     g1->Draw("fcsame");
     g2->Draw("fcsame");

     gStyle->SetOptStat(111);
     gStyle->SetOptFit(1011);
     h->GetPainter()->PaintStat(1011,tot);
     
     h->Draw("e1 same");

     // set the results

     Double_t* result = new Double_t[12];

     for(int i = 0; i<6; i++){
	  result[i]=par[i];
	  result[i+6]=tot->GetParError(i);
     }

     return result;

}