// .L MinimumProblem.C
// MinimumProblem()
//tested on root 5.26/00c
//with roofit3.12

//histograms are yet created with Sumw2() inside

struct FitFunction { 
   FitFunction (TF1 * f) : fFunc(f) {}

   double operator() (const double * x, const double * p) { 
      return p[0]*fFunc->EvalPar(x,&p[1]);
   }
   TF1 * fFunc;
};

RooRealVar CB_mean("CB_mean","CB_mean",120,100,150,"GeV");
RooRealVar CB_dmean("CB_dmean","CB_dmean",0.,-10,10,"GeV");
RooRealVar CB_sigma("CB_sigma","CB_sigma",1.36,0.1,10,"GeV");
RooRealVar CB_alpha("CB_alpha","CB_alpha",1.09,0,2,"GeV" );
RooRealVar CB_n("CB_n","CB_n",6.95,0,200,"");
RooRealVar Gauss_mean("Gauss_mean","Gauss_mean",120,100,150,"GeV");
RooRealVar Gauss_sigma("Gauss_sigma","Gauss_sigma",5.37,1,50,"GeV");
RooRealVar frac_CB("frac_CB","frac_CB",0.97,0,1);

float CB_mean_initial; float error_CB_mean_initial;
float CB_dmean_initial; float error_CB_dmean_initial;
float CB_sigma_initial; float error_CB_sigma_initial;
float CB_alpha_initial; float error_CB_alpha_initial;
float CB_n_initial; float error_CB_n_initial;
float Gauss_mean_initial; float error_Gauss_mean_initial;
float Gauss_sigma_initial; float error_Gauss_sigma_initial;
float frac_CB_initial; float error_frac_CB_initial;

int MinimumProblem()
{
  //the histograms that are opened have yet the option Sumw2()

  void ComputeInvariantMass(TH1F *hist);
  
  using namespace RooFit;
  
  CB_dmean.setConstant(1);
  SaveContext();

  //####################################################################
  RestoreContext();
  TFile *myfile=new TFile("myfile.root");
  myfile->cd();
  ComputeInvariantMass(hist_InvariantMass_TYPE_ggH_mH120_McAtNlo_Jimmy_PileUp2);

  RestoreContext();
  //---
  cout << "##########################" << endl;
  cout << "##########################" << endl;
  cout << "##########################" << endl;
  cout << "##########################" << endl;
  
  TH1F *hist_Total_InvariantMass_H120_PileUp2=new TH1F("hist_Total_InvariantMass_H120_PileUp2","hist_Total_InvariantMass_H120_PileUp2",100,100,150);
  hist_Total_InvariantMass_H120_PileUp2->Sumw2();
  
  //with the following line, the parameters of the fit is the same as for the initial histogram
  //  hist_Total_InvariantMass_H120_PileUp2->Add(hist_InvariantMass_TYPE_ggH_mH120_McAtNlo_Jimmy_PileUp2,1);
  //with the following line, the parameters of the fit is different to the inital histogram
  //  hist_Total_InvariantMass_H120_PileUp2->Add(hist_InvariantMass_TYPE_ggH_mH120_McAtNlo_Jimmy_PileUp2,2);
  hist_Total_InvariantMass_H120_PileUp2->Add(hist_InvariantMass_TYPE_ggH_mH120_McAtNlo_Jimmy_PileUp2,10000);
  
  cout << "bin content(40) : " << hist_InvariantMass_TYPE_ggH_mH120_McAtNlo_Jimmy_PileUp2->GetBinContent(40) << endl;
  cout << "bin error(40) : " << hist_InvariantMass_TYPE_ggH_mH120_McAtNlo_Jimmy_PileUp2->GetBinError(40) << endl;
  
  cout << "bin content(40) : " << hist_Total_InvariantMass_H120_PileUp2->GetBinContent(40) << endl;
  cout << "bin error(40) : " << hist_Total_InvariantMass_H120_PileUp2->GetBinError(40) << endl;
  
  ComputeInvariantMass(hist_Total_InvariantMass_H120_PileUp2);  
}

//##############################################
void ComputeInvariantMass(TH1F *hist) {
  
  hist->Sumw2();
  
  TString myfile_canvas="figures/";
  myfile_canvas+=hist->GetTitle();
  myfile_canvas+=".eps";

  cout << "entries= " << hist->GetEntries() << endl;

  //--------
  TCanvas *canvasRooFit_InvariantMass=new TCanvas("canvasRooFit_InvariantMass","canvasRooFit_InvariantMass",800,600);
  canvasRooFit_InvariantMass->SetRightMargin(0.03);
  canvasRooFit_InvariantMass->SetTopMargin(0.1);
  
    RooRealVar RooRealVar_InvariantMass("RooRealVar_InvariantMass","RooRealVar_InvariantMass",100,150,"GeV"); // TH1 => RooDataSet
  //RooRealVar RooRealVar_InvariantMass("RooRealVar_InvariantMass","RooRealVar_InvariantMass",115,125,"GeV"); // TH1 => RooDataSet
  RooArgSet *observables=new RooArgSet(RooRealVar_InvariantMass);
  //  Bool_t impDens=kFALSE;
  RooDataHist Datahist("Datahist","Datahist",RooRealVar_InvariantMass,RooFit::Import(*hist));//,RooFit::Weight(1./10000));
  //RooDataHist Datahist("Datahist","Datahist",RooRealVar_InvariantMass);
  
  //  RooArgList myavrglist;
  //  myarglist.add(RooRealVar_InvariantMass);
  //  Datahist.importTH1(myarglist,*hist,1.);
  //  importTH1(const RooArgList& vars, TH1& histo, Double_t initWgt, Bool_t doDensityCorrection)

  //i tried with the following lines to see if it can help, but seems not
  Datahist.weightError(RooAbsData::SumW2);
  //Datahist.weightError(RooAbsData::Poisson);
  
  cout << "nb entries of the file imported : " << Datahist.numEntries() << std::endl;
  //=======================
  CB_mean.setVal(120);
  CB_mean.setRange(100,140);
  
  Gauss_mean.setVal(120);
  Gauss_mean.setRange(100,150);
  
  RooFormulaVar* CB_meanFormula=new RooFormulaVar("CB_meanFormula","CB_meanFormula","@0+@1", RooArgList(CB_mean,CB_dmean));
  RooAbsPdf* CB_Pdf=new RooCBShape("CB_Pdf","CB_Pdf",RooRealVar_InvariantMass,*CB_meanFormula,CB_sigma,CB_alpha,CB_n);
  RooAbsPdf* Gauss_Pdf=new RooGaussian("Gauss_Pdf","Gauss_Pdf",RooRealVar_InvariantMass,Gauss_mean,Gauss_sigma);
  RooAddPdf* myPdf_mgg=new RooAddPdf("myPdf_mgg","Signal PDF for H->gamma gamma",RooArgList(*CB_Pdf,*Gauss_Pdf),RooArgList(frac_CB));

  
  cout << "########################################" << endl;
  cout << "########################################" << endl;

  // remove dmean which is constant
  TF1 *f=myPdf_mgg->asTF(RooArgList(RooRealVar_InvariantMass),RooArgList(CB_mean,CB_sigma,CB_alpha,CB_n,Gauss_mean,Gauss_sigma,frac_CB));
  // create the un-normalized function
  
  FitFunction * fitfunc = new FitFunction(f);
  TF1 * f2 = new TF1("fitFunc",fitfunc,100,150,f->GetNpar()+1,"FitFunction");
  f2->SetParName(0,"A");
  f2->SetParameter(0,hist->GetMaximum());
  for (int i = 0; i < f->GetNpar(); ++i) {
     f2->SetParameter(i+1,f->GetParameter(i) );
     f2->SetParName(i+1,f->GetParName(i) );
  }
  f2->SetParLimits(f2->GetParNumber("CB_mean"),CB_mean.getMin(),CB_mean.getMax());
  f2->SetParLimits(f2->GetParNumber("CB_sigma"),CB_sigma.getMin(),CB_sigma.getMax());
  f2->SetParLimits(f2->GetParNumber("CB_alpha"),CB_alpha.getMin(),CB_alpha.getMax());
  f2->SetParLimits(f2->GetParNumber("CB_n"),CB_n.getMin(),CB_n.getMax());
  f2->SetParLimits(f2->GetParNumber("Gauss_mean"),Gauss_mean.getMin(),Gauss_mean.getMax());
  f2->SetParLimits(f2->GetParNumber("Gauss_sigma"),Gauss_sigma.getMin(),Gauss_sigma.getMax());
  f2->SetParLimits(f2->GetParNumber("frac_CB"),frac_CB.getMin(),frac_CB.getMax());

  //f->FixParameter(f.GetParNumber("CB_dmean"),CB_dmean.getVal());
  //TF1 *f=myPdf_mgg->asTF(RooArgList(RooRealVar_InvariantMass));

  // use option 0 to avoid plotting after fit
  hist->Fit(f2,"0");

  //RooFitResult *roofitresult=myPdf_mgg->fitTo(Datahist,RooFit::DataError(RooAbsData::SumW2));
  //i tried with the following to see if if can help but it seems not
  //RooFitResult *roofitresult=myPdf_mgg->fitTo(Datahist,RooFit::SumW2Error(0));
  //RooFitResult *roofitresult=myPdf_mgg->fitTo(Datahist,RooFit::SumW2Error(1));
  
  

  //RooFitResult *roofitresult=myPdf_mgg->chi2FitTo(Datahist,Range(115.,125.));
  
  RooPlot *rooplot_InvariantMass=RooRealVar_InvariantMass.frame(80);
  
  Datahist.plotOn(rooplot_InvariantMass,RooFit::DataError(RooAbsData::SumW2));
  myPdf_mgg->plotOn(rooplot_InvariantMass,LineColor(kBlack),LineWidth(1));
  rooplot_InvariantMass->Draw();
  
  canvasRooFit_InvariantMass->SaveAs(myfile_canvas);

 
  delete canvasRooFit_InvariantMass;
}

void SaveContext() {
  CB_mean_initial=CB_mean.getVal(); error_CB_mean_initial=CB_mean.getError();
  CB_dmean_initial=CB_dmean.getVal(); error_CB_dmean_initial=CB_dmean.getError();
  CB_sigma_initial=CB_sigma.getVal(); error_CB_sigma_initial=CB_sigma.getError();
  CB_alpha_initial=CB_alpha.getVal(); error_CB_alpha_initial=CB_alpha.getError();
  CB_n_initial=CB_n.getVal(); error_CB_n_initial=CB_n.getError();
  Gauss_mean_initial=Gauss_mean.getVal(); error_Gauss_mean_initial=Gauss_mean.getError();
  Gauss_sigma_initial=Gauss_sigma.getVal(); error_Gauss_sigma_initial=Gauss_sigma.getError();
  frac_CB_initial=frac_CB.getVal(); error_frac_CB_initial=frac_CB.getError();
}

void RestoreContext() {
  CB_mean.setVal(CB_mean_initial); CB_mean.setError(error_CB_mean_initial);
  CB_dmean.setVal(CB_dmean_initial); CB_dmean.setError(error_CB_dmean_initial);
  CB_sigma.setVal(CB_sigma_initial); CB_sigma.setError(error_CB_sigma_initial);
  CB_alpha.setVal(CB_alpha_initial); CB_alpha.setError(error_CB_alpha_initial);
  CB_n.setVal(CB_n_initial); CB_n.setError(error_CB_n_initial);
  Gauss_mean.setVal(Gauss_mean_initial); Gauss_mean.setError(error_Gauss_mean_initial);
  Gauss_sigma.setVal(Gauss_sigma_initial); Gauss_sigma.setError(error_Gauss_sigma_initial);
  frac_CB.setVal(frac_CB_initial); frac_CB.setError(error_frac_CB_initial);
}