#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() {
  
  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");

  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);
  RooGaussian gauss("gauss", "gauss", x, mean, sg);
 
  RooRealVar mean1("mean1", "mean1", 0, -5, 5);
  RooRealVar sbw("sbw", "sbw", 2, 0, 10); 
  RooBreitWigner bw("bw", "bw", x, mean1, sbw);
 
  RooFFTConvPdf signal("signal","landau (X) gauss", x, gauss, bw) ;

  signal.fitTo(dh, Range(483.,503.));
  signal.plotOn(frame,DrawOption("L"));
  frame->Draw();

}