#include "RooRealVar.h"
#include "RooCategory.h"
#include "RooFormulaVar.h"
#include "RooPolyVar.h"
#include "RooGaussModel.h"
#include "RooBMixDecay.h"
#include "RooEffProd.h"
#include "RooDataSet.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "RooPlot.h"


using namespace std;
using namespace RooFit;



int Dmd_sens(){
  
  RooRealVar *Bd_propertime = new RooRealVar("Bd_propertime_corrected","t",0.001,12, "ps");//_corrected
  RooRealVar *TAU_Bd        = new RooRealVar("#tau_{Bd}","#tau_{Bd}",1.525,0.,2. , "ps");
  RooRealVar *dmd           = new RooRealVar("#Delta m_{d}","##Delta m_{d}",0.507,0.,1., "ps");
  RooRealVar *mistag        = new RooRealVar("mistag","mistag",0.365 , 0. , 0.5);
  RooRealVar *Delta_mistag  = new RooRealVar("Delta_mistag","Delta_mistag",0.);//,-1.,1.);
 
  RooCategory *TagFlav = new RooCategory("TagFlav_true","B0/B0bar mixing state") ;
  TagFlav->defineType("B0",-1) ;
  TagFlav->defineType("B0bar",1) ;
  TagFlav->defineType("untagged",0) ;
    
  RooCategory *mixState = new RooCategory("mixState_reco","B0/B0bar mixing state") ;
  mixState->defineType("mixed",-1);
  mixState->defineType("unmixed",1);
  mixState->defineType("untagged",0);

  RooRealVar *alpha = new RooRealVar("#alpha","#alpha", 2.2 , 0. ,5. );
  RooRealVar *beta  = new RooRealVar("#beta" ,"#beta" , 2.4 , 0. ,5. );
  beta->setConstant(kTRUE);
  alpha->setConstant(kTRUE);
  
  RooFormulaVar *Acc_Function
    = new RooFormulaVar("Acc_Function","Acc_Function",
                        "(@0>0)*pow(@0*@1,@2)/(1+pow(@0*@1,@2))", // *(1+@3*@0)
                        RooArgList(*Bd_propertime,*alpha,*beta));// ,*slope
  
  RooRealVar *p0_mean = new RooRealVar("p0_mean","p0_mean"   , 0.021 , -1 , 1 ) ;
  RooRealVar *p1_mean = new RooRealVar("p1_mean","p1_mean"   ,-0.047 , -1 , 1 ) ;
  RooRealVar *p0_sigma = new RooRealVar("p0_sigma","p0_sigma", 0.076  ,  0 , 2 ) ;
  RooRealVar *p1_sigma = new RooRealVar("p1_sigma","p1_sigma", 0.148  ,  0 , 2 ) ;
  
  p0_mean->setConstant(kTRUE);
  p1_mean->setConstant(kTRUE);
  p0_sigma->setConstant(kTRUE);
  p1_sigma->setConstant(kTRUE);
 
  RooPolyVar mean("mean","mean",*Bd_propertime,RooArgSet(*p0_mean,*p1_mean));
  RooPolyVar sigma("sigma","sigma",*Bd_propertime,RooArgSet(*p0_sigma,*p1_sigma));
  RooGaussModel res("res","res",*Bd_propertime,mean,sigma);
    
  RooBMixDecay *Bd2DstarMuNu_Theopdf
    = new RooBMixDecay("Bd2DstarMuNu_Theopdf","Bd2DstarMuNu_Theopdf",
                       *Bd_propertime,*mixState,*TagFlav,*TAU_Bd,*dmd, // Physics parameters
                       *mistag,*Delta_mistag,                          // Mistag parameters
		       res,
                       RooBMixDecay::SingleSided);
  
  
  RooEffProd *Bd2DstarMuNu_pdf_acc =
    new RooEffProd("Bd2DstarMuNu_pdf_acc", "Bd2DstarMuNu_pdf_acc",
		   *Bd2DstarMuNu_Theopdf, 
		   *Acc_Function);
  
  RooEffProd *Total_pdf_acc
    = new RooEffProd("Total_pdf_acc","Total_pdf_acc",
		     *Bd2DstarMuNu_Theopdf,*Acc_Function);
  
  
  RooDataSet *data = Total_pdf_acc->generate(RooArgSet( *Bd_propertime , *mixState ,*TagFlav ) , 200000); 
  RooFitResult *result = Total_pdf_acc->fitTo(*data,NumCPU(24),PrintLevel(-1),Save(true)); // 
  result->Print("v");
  
  TCanvas *c_propertimefit= new TCanvas("c_propertimefit","c_propertimefit");
  RooPlot* frame_t_all = Bd_propertime->frame(Title("All Events")) ;
  data->plotOn(frame_t_all);
  Total_pdf_acc->plotOn(frame_t_all);
  Total_pdf_acc->paramOn(frame_t_all);
  frame_t_all->Draw();
 
  return 0;
}