#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooLandau.h"
#include "RooFFTConvPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooChi2Var.h"
#include <RooCmdArg.h>

using namespace RooFit;

#include <cmath>

void RooFit_MASS() {

#if 0
  TChain chain("h1");
  chain.Add("SelppM_new_0_9.root");
  chain.Add("SelppM_new_10_29.root");
  chain.Add("SelppM_new_30_50.root");

  TH1D *hma1 = new TH1D("hma1","M_invariant for 2-body decay", 200, 475., 512.);

  TCut cutcos="sqrt(pkvt[0][0]*pkvt[0][0]+pkvt[0][1]*pkvt[0][1])/sqrt(pkvt[0][0]*pkvt[0][0]+pkvt[0][1]*pkvt[0][1]+pkvt[0][2]*pkvt[0][2])>=0.9&&sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1])/sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1]+ppvt[0][2]*ppvt[0][2])>=0.9";    
  TCut spion="pstarpvt[0]<220.";
  TCut fkin = "sqrt(rkvt[0][0]*rkvt[0][0]+rkvt[0][1]*rkvt[0][1])>30"; 
  TCut cutsign = "chk[0]*cht[0]*flip[0]==1";
  TCut add = "sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1])>= 140.&&sqrt(ppvt[0][0]*ppvt[0][0]+ppvt[0][1]*ppvt[0][1])<=240.&&sqrt(pkvt[0][0]*pkvt[0][0]+pkvt[0][1]*pkvt[0][1])>80.&&Minvariant[0]<=520.&&Minvariant[0]>=470.";


  chain.Draw("Minvariant[0]>>hma1", cutcos&&spion&&fkin&&cutsign&&add&&"nt==1&&nkg==1");
#endif

    auto f = TFile::Open("histogram.root");
  auto c1 = (TCanvas*) f->Get("c1");
  auto hma1 = (TH1*) c1->FindObject("hma1")->Clone();
  //f->Close();

  RooRealVar x("x", "x", 475., 512.);
  RooDataHist dh("dh", "dh", x, Import(*hma1));

  RooPlot *frame = x.frame(Title("M_invariant for 2-body decay"));
  dh.plotOn(frame);

  RooRealVar mean("mean", "mean", 493.677, 490, 500);
  RooRealVar sg("sg", "sg", 1, 0, 10);
  RooRealVar mean1("mean1", "mean1", 0);
  RooGaussian gauss("gauss", "gauss", x, mean1, sg);
 
  RooRealVar sbw("sbw", "sbw", 2, 0, 10); 
  RooBreitWigner bw("bw", "bw", x, mean, sbw);
 
  RooFFTConvPdf signal("signal","landau (X) gauss", x, bw, gauss) ;

  x.setBins(100000, "cache");
  //signal.setBufferFraction(0.5);


  // RooVoigtian signal("signal","voigt",x, mean, sbw, sg);

  auto r = signal.fitTo(dh, Range(483.,503.), Save(1), Minimizer("Minuit2") /*, Minos(RooArgSet(mean)) */ );
  signal.plotOn(frame,DrawOption("L"));
  new TCanvas(); 
  frame->Draw();
  r->Print();
}