{

  gSystem->Load("libRooFit");
  using namespace RooFit;

  TChain *ch = new TChain("tree");
  ch->Add("/home/deepak/BsToPhiMuMu/DATA_DoubleMuon/cut1/sel_BsToPhiMuMu_Double_2016*.root");
  ch->Add("/home/deepak/BsToPhiMuMu/DATA_Charmonium/Rereco/cut1/sel_BsToPhiMuMu_*.root");

  TTree *tr = ch;

  int nentries_ = tr->GetEntries();
  cout << "\n=> total entries in the data BsToPhiMuMu tree = " << nentries_ << endl;


  RooRealVar Bmass("Bmass","M(B_{s}) [GeV/c^{2}]", 4.7, 6.0);
  RooRealVar Mumumass("Mumumass","M^{#mu#mu}", 0., 10.);
  RooRealVar Mumumasserr("Mumumasserr","Error of M^{#mu#mu}", 0., 10.);
  RooRealVar Triggers("Triggers","HLT", 0, 1);
  RooRealVar Phimass("Phimass","M^{#phi}", 0., 2.);
  
  RooArgSet s(Bmass, Mumumass,Mumumasserr,Triggers,Phimass);

  RooDataSet data("data","dataset with Bmass", ch, RooArgSet(s));

  TCut trig   = "Triggers==1";
  TCut phim = "Phimass > 1.01 && Phimass < 1.03";

  TCut jpR = "(Mumumass < 3.096916-5.5*Mumumasserr || Mumumass > 3.096916+3.5*Mumumasserr)";
  TCut psR = "(Mumumass < 3.686109-3.5*Mumumasserr || Mumumass > 3.686109+3.5*Mumumasserr)";

  //optimized anti-radiation cuts
  TCut AR = "abs(Bmass - Mumumass - 2.270) > 0.170 && abs(Bmass - Mumumass - 1.681) > 0.080";

  RooDataSet *redData = (RooDataSet*)data.reduce(trig && phim  && jpR && psR && AR);

  RooRealVar  mean("mean","mean of the gaussians", 5.367, 5.1, 5.6);
  RooRealVar  sigma1("sigma1","sigma1 of the gaussian1",0.03206);//fixed sigma
  RooGaussian sig1("sig1","Signal component 1",Bmass,mean,sigma1) ;

  RooRealVar  sigma2("sigma2","sigma2 of the gaussian2",0.11837);//fixed sigma
  RooGaussian sig2("sig2","Signal component 2",Bmass,mean,sigma2) ;

  RooRealVar sigfrac("sigfrac","fraction of component  in signal",0.88208) ;//fixed signal fraction
  RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sigfrac) ;

  //  RooRealVar lambda("lambda","slope",-5.0,3.5);
  //RooExponential bkgExpo("bkgExpo","exponential PDF", Bmass, lambda);

  RooRealVar poly_c1("poly_c1","coefficient of linear polynomial",0, -5, 5);
  RooPolynomial linpoly("linpoly","linear polynomial",Bmass,poly_c1);

  RooRealVar poly_c2("poly_c2","coefficient of polynomial",0, -5, 5);
  RooPolynomial poly("poly","polynomial for datasideband",Bmass,poly_c2);
  
  
  RooRealVar nsigG("sigG","signal events", 0, 2459);
  RooRealVar npolyd("npolyd","number of bkg", 0, 2459);
  RooRealVar npoly("npoly","number of bkg events", 0, 2459);

  RooAddPdf model("model","",RooArgSet(sig,linpoly),RooArgSet(nsigG,npoly));
  model.Print("t");

  Bmass.setRange("signal",5.1,5.6);
  RooFitResult* fitres = model.fitTo(*redData, Save(true),Range("signal") );  

  Bmass.setRange("back",4.7,5.6);
  RooFitResult* pol = poly.fitTo(*redData, Save(true),Range("back") );
  
  RooPlot *xframe = Bmass.frame(Bins(65) );
  xframe->SetTitle("");
  redData->plotOn(xframe);
  model.plotOn(xframe);

  /*  model.plotOn(xframe, Components("sig"),LineStyle(kDashed), LineColor(kRed),FillColor(kGreen), FillStyle(3000), DrawOption("F") );
  model.plotOn(xframe, Components("bkgExpo"), LineStyle(kDotted),LineColor(kRed),FillColor(kBlue), LineColor(kRed) , FillStyle(3000),DrawOption("F"));
  model.plotOn(xframe, Components("linpoly"), FillColor(kYellow), LineColor(kYellow) , FillStyle(3000),DrawOption("F"));
  */

  model.plotOn(xframe, Components("sig"), LineStyle(kDashed), LineColor(kRed), FillColor(kGreen), FillStyle(3000), DrawOption("F") );
  model.plotOn(xframe, Components("linpoly"),FillColor(kRed),FillStyle(3000),DrawOption("F"));
  model.plotOn(xframe, Components("poly"), FillColor(kYellow), LineColor(kYellow) , FillStyle(3000),DrawOption("F")); 
  
  TCanvas *c = new TCanvas("c","c", 600,600);
  gPad->SetLeftMargin(0.15);
  xframe->GetYaxis()->SetTitleOffset(1.7);
  xframe->Draw();

  TCanvas *c = new TCanvas("c","c", 600,600);
  gPad->SetLeftMargin(0.15);
  xframe->GetYaxis()->SetTitleOffset(1.7);
  xframe->Draw();

  TPaveText* paveText = new TPaveText(0.62,0.70,0.80,0.88,"NDC");
  paveText->SetBorderSize(0.0);
  paveText->SetFillColor(kWhite);
  paveText->SetFillStyle(0);
  paveText->SetTextSize(0.02);
  paveText->AddText(Form("nsigG = %.0f #pm %.0f evts", nsigG.getVal(), nsigG.getError()));
  paveText->AddText(Form("npoly = %.0f #pm %.0f evts", npoly.getVal(), npoly.getError()));
  paveText->AddText(Form("npolyd = %.0f #pm %.0f evts", npolyd.getVal(), npolyd.getError()));
  paveText->AddText(Form("mean = %.5f #pm %.5f GeV", mean.getVal(), mean.getError()));
  paveText->AddText(Form("sigma1 = %.5f #pm %.5f GeV", sigma1.getVal(), sigma1.getError()));
  paveText->AddText(Form("sigma2 = %.5f #pm %.5f GeV", sigma2.getVal(), sigma2.getError()));
  paveText->AddText(Form("sigfrac = %.5f #pm %.5f GeV", sigfrac.getVal(), sigfrac.getError()));

  if (fitres != NULL){
    if ((fitres->covQual() == 3) && (fitres->status() == 0)){
      if (paveText != NULL) paveText->AddText("Fit Status: GOOD");
    } else {
      if (paveText != NULL) paveText->AddText("Fit Status: BAD");
    }
  }

  paveText->Draw();

  TLine *l1 = new TLine(5.367-7.0*0.032,0,5.367-7.0*0.032,100);
  l1->SetLineColor(kRed);
  l1->Draw();

  TLine *l2 = new TLine(5.367-4.5*0.032,0,5.367-4.5*0.032,100);
  l2->SetLineColor(kRed);
  l2->SetLineStyle(10);
  l2->Draw();

  TLine *l3 = new TLine(5.367+4.5*0.032,0,5.367+4.5*0.032,100);
  l3->SetLineColor(kBlue);
  l3->Draw();

  TLine *l4 = new TLine(5.367+7.0*0.032,0,5.367+7.0*0.032,100);
  l4->SetLineColor(kBlue);
  l4->Draw();

  
  TLatex *LumiTex1 = new TLatex();
  LumiTex1->SetTextColor(kBlack);
  LumiTex1->SetNDC(true);

  TLatex *LumiTex2 = new TLatex();
  LumiTex2->SetTextColor(kBlack);
  LumiTex2->SetNDC(true);

  TLatex *LumiTex3 = new TLatex();
  LumiTex3->SetTextColor(kBlack);
  LumiTex3->SetNDC(true);

  TString LumiName1 = "CMS" ;
  TString LumiName2 = "#it{Preliminary}" ;
  TString LumiName3 = "35.9fb^{-1} (13TeV)" ;

  LumiTex1->DrawLatex(0.18,0.9,LumiName1);
  LumiTex2->DrawLatex(0.3,0.9,LumiName2);
  LumiTex3->DrawLatex(0.6,0.9,LumiName3);


  c->Update();
  
  c->SaveAs("Bmass_full.png");
  
}