void bgEstTemplateShapeTest()
{
  const unsigned numToyExperiments = 1000;

  TH1* histogramI   = testTemplateShape(10000,   100, "I", numToyExperiments);
  TH1* histogramII  = testTemplateShape(  100, 10000, "II", numToyExperiments);
  TH1* histogramIII = testTemplateShape(  100,   100, "III", numToyExperiments);
  TH1* histogramIV  = testTemplateShape(10000, 10000, "IV", numToyExperiments);

  TCanvas* canvas = new TCanvas("canvas", "canvas", 1, 1, 800, 600);
  canvas->SetFillColor(10);
  canvas->SetBorderSize(2);

  TPostScript* ps = new TPostScript("bgEstTemplateShapeTest.ps", 112);

  TLegend legend(0.53, 0.73, 0.88, 0.89);
  legend.SetBorderSize(0);
  legend.SetFillColor(0);

  double yMax = 0.;
  if ( histogramI->GetMaximum()   > yMax ) yMax = histogramI->GetMaximum();
  if ( histogramII->GetMaximum()  > yMax ) yMax = histogramII->GetMaximum();
  if ( histogramIII->GetMaximum() > yMax ) yMax = histogramIII->GetMaximum();
  if ( histogramIV->GetMaximum()  > yMax ) yMax = histogramIV->GetMaximum();
  
  histogramI->SetMaximum(1.3*yMax);
  histogramI->SetMinimum(0.);
  histogramI->SetStats(false);
  histogramI->SetLineColor(2);
  histogramI->SetLineWidth(2);
  legend.AddEntry(histogramI, "I", "l");
  histogramI->Draw();

  histogramII->SetLineColor(3);
  histogramII->SetLineWidth(2);
  legend.AddEntry(histogramII, "II", "l");
  histogramII->Draw("same");
  
  histogramIII->SetLineColor(4);
  histogramIII->SetLineWidth(2);
  legend.AddEntry(histogramIII, "III", "l");
  histogramIII->Draw("same");

  histogramIV->SetLineColor(7);
  histogramIV->SetLineWidth(2);
  legend.AddEntry(histogramIV, "IV", "l");
  histogramIV->Draw("same");

  legend.Draw();

  canvas->Update();
  
  ps->NewPage();

  delete ps;

  delete canvas;

  delete histogramI;
  delete histogramII;
  delete histogramIII;
  delete histogramIV;
}

TH1* sampleDistribution(const char* name, unsigned numEvents, 
			float coeff, float offset1, float slope1, float offset2, float slope2, 
			bool normalize)
{
  TH1* histogram = new TH1F(name, name, 20, 0., 10.);
  
  static TRandom3 rnd;
  
  for ( unsigned iEvent = 0; iEvent < numEvents; ++iEvent ) {
    bool histogram_filled = false;
    while ( !histogram_filled ) {
      float u_x = 10*rnd.Uniform();
      float u_y = 10*rnd.Uniform();
      if ( u_y <= (coeff*(offset1 + slope1*u_x) + (1 - coeff)*(offset2 + slope2*u_x)) ) {
	histogram->Fill(u_x);
	histogram_filled = true;
      }
    }
  }

  if ( normalize ) histogram->Scale(1./histogram->Integral());

  return histogram;
}

TH1* testTemplateShape(unsigned numEvents1, unsigned numEvents2, const char* label, unsigned numToyExperiments)
{
  std::cout << "testing RooFit with event statistics (1) = " << numEvents1 << ", (2) = " << numEvents2 << "..." << std::endl;

  TString histogramName = TString("FitResults").Append(label);
  TH1* histogram = new TH1F(histogramName.Data(), histogramName.Data(), 202, -0.51, +1.51);

  for ( unsigned iToyExperiment = 0; iToyExperiment < numToyExperiments; ++iToyExperiment ) {
    std::cout << "iToyExperiment = " << iToyExperiment << std::endl;

    RooRealVar x("x", "x", 0., 10.);

    std::cout << " preparing modelTemplate1..." << std::endl;
    RooRealVar offset1("offset1", "offset1", 6., 0., 10.);
    RooRealVar slope1("slope1", "slope1", +0.2, -1., +1.);
    RooGenericPdf modelPdf1("modelPdf1","offset1 + slope1*x", RooArgList(offset1, slope1, x));
    //RooDataHist* data1 = modelPdf1.generateBinned(RooArgSet(x), numEvents1);
    TH1* data1 = sampleDistribution("data1", numEvents1, 1.0, 6., +0.2, 0., 0., true);
    RooDataHist* template1 = new RooDataHist("template1", "template1", RooArgList(x), data1);
    RooHistPdf modelTemplate1("modelTemplate1", "modelTemplate1", RooArgSet(x), *template1);
    
    std::cout << " preparing modelTemplate2..." << std::endl;
    RooRealVar offset2("offset2", "offset2", 8., 0., 10.);
    RooRealVar slope2("slope2", "slope2", -0.2, -1., +1.);
    RooGenericPdf modelPdf2("modelPdf2","offset2 + slope2*x", RooArgList(offset2, slope2, x));
    //RooDataHist* data2 = modelPdf2.generateBinned(RooArgSet(x), numEvents2);
    TH1* data2 = sampleDistribution("data2", numEvents2, 0.0, 0., 0., 8., -0.2, true);
    RooDataHist* template2 = new RooDataHist("template2", "template2", RooArgList(x), data2);
    RooHistPdf modelTemplate2("modelTemplate2", "modelTemplate2", RooArgSet(x), *template2);

    std::cout << " preparing sumData..." << std::endl;
    RooRealVar coeffModel("coeffModel", "coeffModel", 0.75, 0., 1.);
    RooAddPdf sumModel("sumModel", "coeff*modelPdf1 + (1 - coeff)*modelPdf2", modelPdf1, modelPdf2, coeffModel);
    //RooDataHist* sumData = sumModel.generateBinned(RooArgSet(x), 1000);
    TH1* sumData_aux = sampleDistribution("sumData_aux", 1000, 0.75, 6., +0.2, 8., -0.2, false);
    RooDataHist* sumData = new RooDataHist("sumData", "sumData", RooArgList(x), sumData_aux);
    
    std::cout << " executing fit..." << std::endl;
    RooRealVar coeffTemplate("coeffTemplate", "coeffTemplate", 0.50, -0.5, +1.5);
    RooAddPdf sumTemplate("sumTemplate", "coeff*modelTemplate1 + (1 - coeff)*modelTemplate2", modelTemplate1, modelTemplate2, coeffTemplate);
    sumTemplate.fitTo(*sumData);

    std::cout << "--> coeffTemplate = " << coeffTemplate.getVal() << std::endl;
    histogram->Fill(coeffTemplate.getVal());

    std::cout << " deleting temporary objects..." << std::endl;
    delete template1;
    delete data1;

    delete template2;
    delete data2;

    delete sumData;
    delete sumData_aux;
  }

  return histogram;
}