//
// 30-11-09
//  
// closure test: take signal region template and fit it back
//-----------------------------
#include <iomanip>

using namespace RooFit;

void fit_closure_3(double v1 = 1.0, double v2 = 2.0){


  gStyle->SetOptStat(1111111);
  gStyle->SetStatW(0.35);

  double fit_from = 0;
  double fit_upto = 0.2;



  TCanvas c0("c0","fit",3*400,500);
  c0.Divide(3,1);

  const double tiny = 1e-6;

  ///------------------------------------
  ///  Get template from control sample
  ///------------------------------------


  double ylow = 0;///5.8;//0;
  double yup  = 5;//13.2;//4;
  if(v1>9) { ylow = 8; yup=15; }

  // make up a 2-bin histogram
  TH1F *h = new TH1F("h","RelIso",2,0,0.2);
  //double v1 = 1;//9.9;//8.93624;//1;
  //double v2 = 2.;//9.9;// ;9.58708;//3;
  h->SetBinContent(1,v1);
  h->SetBinContent(2,v2);

  c0.cd(1);
  gPad->SetTicky();
  h->SetLineColor(kBlack);
  h->GetYaxis()->SetRangeUser( ylow, yup );
  h->GetXaxis()->SetNdivisions( 102 );
  h->SetTitle("Mock data histogram");
  h->SetMarkerSize(2.5);
  h->DrawCopy("h text0");




  double ntrue = h->Integral(1,2);
  cout << "sum of bin 1 and 2 = " << ntrue << endl;



  // create HistPdf
  RooRealVar reliso("reliso","reliso", 0, 0.2);
  reliso.setBins(2);
  RooDataHist rh_qcd("rh_qcd","qcd", reliso, h);


  RooHistPdf pdf_qcd( "pdf_qcd", "QCD pdf ", reliso, rh_qcd, 0);

  // create model
  //  double  nlimit  = 2*ntrue;
  double  nlimit  = 10;
  if( v1>=10 ) nlimit = 100;
  RooRealVar nqcd("nfit", "number of events", 0, nlimit, "event");
  RooAddPdf model("model","model", pdf_qcd, nqcd);

 
  ///--------------------------
  ///  Fit signal sample: "all"
  ///---------------------------

  // Make a copy of the original histogram to fit
  //  TH1F *all = (TH1F*)h->Clone();

  // Prepare data in *fit range*
  //  all->GetXaxis()->SetRangeUser(fit_from,fit_upto);
  RooDataHist data("data","data in signal sample",reliso, h);

  // By default is extended ML.
  model.fitTo( data, SumW2Error(kFALSE) );// ROOT >= 5.24
   //model.fitTo( data );// ROOT 5.22


  /*
[#0] WARNING:InputArguments -- RooAbsPdf::fitTo(model) WARNING: a likelihood fit is request of what appears to be weighted data. 
       While the estimated values of the parameters will always be calculated taking the weights into account, 
       there are multiple ways to estimate the errors on these parameter values. You are advised to make an 
       explicit choice on the error calculation: 
           - Either provide SumW2Error(kTRUE), to calculate a sum-of-weights corrected HESSE error matrix 
             (error will be proportional to the number of events)
           - Or provide SumW2Error(kFALSE), to return errors from original HESSE error matrix
             (which will be proportional to the sum of the weights)
       If you want the errors to reflect the information contained in the provided dataset, choose kTRUE. 
       If you want the errors to reflect the precision you would be able to obtain with an unweighted dataset 
       with 'sum-of-weights' events, choose kFALSE.
  */




  cout << "\ntrue number = " << ntrue << endl;
  cout << "fit resutls" << endl;
  nqcd.Print();
  double nfit = nqcd.getVal();
  cout << endl;


  // plot with DataError: Poisson (Def) or SumW2
  c0.cd(2);


  h->SetStats(0);
  h->SetTitle("Binned ML Fit");
  h->GetXaxis()->SetNdivisions( 102 );
  h->DrawCopy();
  gPad->SetTicky();

  RooPlot *isoF1 = reliso.frame();
  int m1 = 20;//22; //triangle up
  int m2 = 21;//23; //triangle down
  data.plotOn(isoF1, DataError(RooAbsData::SumW2), 
	      LineWidth(1),
	      LineStyle(kDashed), 
	      LineColor(kBlue), 
	      MarkerColor(kBlue),
	      MarkerStyle(m1));
  data.plotOn(isoF1, DataError(RooAbsData::Poisson),
	      LineWidth(1),
	      LineStyle(kDotted), 
	      LineColor(kGreen+2), 
	      MarkerColor(kGreen+2),
	      MarkerStyle(m2) );

  model.plotOn(isoF1, 
	       LineWidth(1),
	       LineColor(kRed), 
	       LineStyle(kDashDotted)
	       );

  model.paramOn(isoF1, //Label("fit result"),
		Layout(0.4,0.85,0.93));
  //  isoF1->SetTitle("DataError=Poisson");
  isoF1->SetYTitle("");
  //isoF1->SetAxisRange(ylow,yup,"y");
  isoF1->GetXaxis()->SetNdivisions(102);
  isoF1->Draw("same");



  c0.cd(3);
  // 1. original histo
  h->SetStats(0);
  h->SetTitle("Binned ML fit, model scaled to nfit");
  //  h->GetXaxis()->SetNdivisions( 102 );
  h->DrawCopy();
  gPad->SetTicky();


  // RooDataHist
  RooPlot *isoF2 = reliso.frame();
  // 2a) original (SumW2)
  // 2b) rounded (Poisson)
  data.plotOn(isoF2, DataError(RooAbsData::SumW2), 
	      LineWidth(1),
	      LineStyle(kDashed), 
	      LineColor(kBlue), 
	      MarkerColor(kBlue),
	      MarkerStyle(m1) );
  data.plotOn(isoF2, DataError(RooAbsData::Poisson), 
	      LineWidth(1),
	      LineStyle(kDotted), 
	      LineColor(kGreen+2), 
	      MarkerColor(kGreen+2),
	      MarkerStyle(m2) );
  
  model.plotOn(isoF2, 
	       LineWidth(1), 
	       LineColor(kRed), 
	       LineStyle(kDashDotted), 
	       Normalization(nfit, RooAbsReal::NumEvent)
	       );

  model.paramOn(isoF2, //Label("fit result"),
		Layout(0.4,0.85,0.93));
  //  isoF2->SetTitle("DataError=Poisson");
  // isoF2->SetAxisRange(ylow,yup,"y");
  isoF2->GetXaxis()->SetNdivisions(102);
  isoF2->SetYTitle("");
  isoF2->Draw("same");



  //  TH1F *rh_qcd2 = (TH1F*)rh_qcd.createHistogram("rh_qcd2",reliso);
  TLegend leg(0.2,0.65,0.9,0.85);
  leg.AddEntry(h,"mock data histogram","l");

  TH1F *t2 = new TH1F("t2","t2",1,0,1);  
  TH1F *t3 = new TH1F("t3","t3",1,0,1);
  TH1F *t4 = new TH1F("t4","t4",1,0,1);
  t2->SetLineColor(kBlue);
  t2->SetMarkerColor(kBlue);
  t2->SetLineStyle(kDashed);
  t2->SetMarkerStyle(m1);

  t3->SetLineColor(kGreen+2);
  t3->SetMarkerColor(kGreen+2);
  t3->SetLineStyle(kDotted);
  t3->SetMarkerStyle(m2);

  t4->SetLineColor(kRed);
  t4->SetMarkerColor(kRed);
  t4->SetLineStyle(kDashDotted);

  leg.AddEntry(t2,"RooDataHist SumW2 (orig)","lp");
  leg.AddEntry(t3,"RooDataHist Poisson (rounded)","lp");

  TLegend leg2(leg);
  leg2.AddEntry(t4,"model (norm to nfit)","l ");

  leg.AddEntry(t4,"model (norm to 'rounded')","l ");

  c0.cd(2); leg.Draw("same");
  c0.cd(3); leg2.Draw("same");

  bool oneDecimal_v1 = int(v1*100)%10 == 0; //19.9 - 19 =  0.9 > 0.1
  bool twoDecimal_v1 = int(v1*100)%10 > 0; //1.99 -> 19.9 -> 0.9 > 0.1

  bool oneDecimal_v2 = int(v2*100)%10 == 0; //19.9 - 19 =  0.9 > 0.1
  bool twoDecimal_v2 = int(v2*100)%10 > 0; //1.99 -> 19.9 -> 0.9 > 0.1

  //cout << "oneDecimal = " << oneDecimal_v1  << endl;
  //cout << "twoDecimal = " << twoDecimal_v1  << endl;

  char *rootver = "5.25.04";
  if(oneDecimal_v1 && oneDecimal_v2) //eg 1,2
    c0.SaveAs( Form("root_%s/closure_3__%.1f_%.1f.gif",rootver,v1,v2) );

  else if(twoDecimal_v1 && twoDecimal_v2)  //eg 1.99,2.99
    c0.SaveAs( Form("root_%s/closure_3__%.2f_%.2f.gif",rootver,v1,v2) );

  else if(oneDecimal_v1 && twoDecimal_v2) //eg 1.0, 1.99
    c0.SaveAs( Form("root_%s/closure_3__%.1f_%.2f.gif",rootver,v1,v2) );

  else if(twoDecimal_v1 && oneDecimal_v2) //eg 1.99,1.0
    c0.SaveAs( Form("root_%s/closure_3__%.2f_%.1f.gif",rootver,v1,v2) );
    
}