//////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH2.h"
#include "RooFitResult.h"
#include "TStyle.h"
#include "TDirectory.h"
#include "RooGExpModel.h"
#include "Belle2Style.h"
using namespace std;
using namespace RooFit;
void Resolution_function_perevent_post(){
TChain* ch = new TChain("h1_nomva");
 ch->Add("signalmc2.root");
 SetBelle2Style();
 // Fit observables
  RooRealVar res("res", "res",-10,10);
  RooRealVar dterr("dterr", "per-event error on dt", 0.0, 10);
  //Filling the datset
  RooDataSet* data = new RooDataSet("data", "data",RooArgSet(res));  // no of entries/events in the datafile
  Int_t nevt = (int)ch->GetEntries();
  double  DeltaT,DeltaTErr,mcDeltaT,Mbc_mod,deltaE,mySig,BDT,qr;
  ch->SetBranchAddress("DeltaT",&DeltaT);
  ch->SetBranchAddress("DeltaTErr",&DeltaTErr);
  ch->SetBranchAddress("mcDeltaT",&mcDeltaT);
  ch->SetBranchAddress("mod_mbc",&Mbc_mod);
  ch->SetBranchAddress("deltaE",&deltaE);
  ch->SetBranchAddress("isSignal",&mySig);
  ch->SetBranchAddress("BDT",&BDT);
  ch->SetBranchAddress("qrOutput__boFBDT__bc",&qr);
// Import the data points
for(int i=0; i<nevt;i++)
    {

      ch->GetEntry(i);

     if (mySig==1 and deltaE>-0.3 and deltaE<0.3 and Mbc_mod>5.24 and Mbc_mod<5.289 and BDT>0.91 and DeltaTErr<2.5 and abs(qr)>0.1 and abs(qr)<0.25){

              double residual=DeltaT-mcDeltaT;
              if(residual>-10 and residual <10){
              res.setVal(residual);
	      dterr.setVal(DeltaTErr);
            // just for checking things are ok or not
             data->add(RooArgSet(res));
	      }

     }

    }
// For plotting  Deltat Resolution
  RooPlot *res_frame = data->plotOn(res.frame(40));
   //variables and parameters
   //--------------------Resolution Modeling----------------------------------//
   //------  Gaussian + RooGExpModel_normal + RooGExpModel_flipped---------//
   // Main Gaussian mean
   RooRealVar mean_main("mean_main","#mu_{main}", -1.16849e-01, -1,1);
   // Tail Gaussian mean 
   RooRealVar mean_tail("mean_tail","#mu_{tail}", 8.49497e-01, -1,1);
   RooRealVar rlifeL("rlifeL","C_{L}",1, 0.2, 3);
   RooRealVar rlifeR("rlifeR","C_{R}",1, 0.2, 3);
   RooRealVar fmain("fmain", "fmain", 3.66203e-01, 0.1, 1);
   RooRealVar ftail("ftail", "ftail", 5.48298e-01, 0.1, 1);

   ////////---------CORE PART---------------------////////////////
   // Main Gausiian function
   RooGaussModel core("core","core", res, mean_main,dterr);
   //RooGExpModel_normal
  RooGExpModel normal ("normal","normal",res, mean_tail,dterr, rlifeR, RooConst(1.), RooConst(1.), RooConst(1.), kTRUE,RooGExpModel::Normal);
 ////RooGExpModel_flipped
  RooGExpModel flipped ("flipped","flipped",res, mean_tail,dterr, rlifeL, RooConst(1.), RooConst(1.), RooConst(1.), kTRUE,RooGExpModel::Flipped);
  RooAddPdf resolution("resolution","resolution",RooArgList(core,normal,flipped), RooArgList(fmain,ftail));
//////////////////Resoloution Model complete////////

//fitining to the data points
  resolution.fitTo(*data,Save(kTRUE),ConditionalObservables(dterr));
  ///// Style pull ploting
  TCanvas * a1 = new TCanvas("a1","a1",600,600);
  Double_t xlow, ylow, xup, yup;
   a1->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
   a1->Divide(1,2);
   TVirtualPad *upPad = a1->GetPad(1);
   upPad->SetPad(xlow,ylow+0.30*(yup-ylow),xup,yup);

   TVirtualPad *dwPad = a1->GetPad(2);
   dwPad->SetPad(xlow,ylow,xup,ylow+0.35*(yup-ylow));
   dwPad->SetBottomMargin(0.3);
   a1->Update();
   a1->cd(1);
   resolution.plotOn(res_frame);
  resolution.paramOn(res_frame,Format("NEU", AutoPrecision(1)));
   cout << "chi^2_res = " << res_frame->chiSquare() << endl ;
  //pull calculation
   RooPlot* res_zframe = res.frame();
   RooHist* res_hpull = res_frame->pullHist();
   res_hpull->SetFillStyle(1001);
   res_hpull->SetFillColor(1);
    for(int i=0;i<res_hpull->GetN();++i) res_hpull->SetPointError(i,0.0,0.0,0.0,0.0);
   res_frame->Draw();
   //pull plot
  res_zframe->addPlotable(res_hpull,"B");
  a1->cd(2);
  res_zframe->Draw();
}