{

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

  TChain *ch = new TChain("tree");
  ch->Add("/home/deepak/DataFiles_20_04_2018/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}]", 5.1, 5.6);
  RooRealVar Bmass("Bmass","M(B_{s}) [GeV/c^{2}]", 5.2, 5.6);
  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 bm = "Bmass > 5.1 && Bmass < 5.6";
  TCut bm = "Bmass > 5.2 && Bmass < 5.6";
  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)";

  TCut AR = "abs(Bmass - Mumumass - 2.270) > 0.170 && abs(Bmass - Mumumass - 1.681) > 0.080";

  RooDataSet *redData = (RooDataSet*)data.reduce(trig && bm && 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
  //RooRealVar  sigma1("sigma1","sigma1 of the gaussian1",0,1);
  RooGaussian sig1("sig1","Signal component 1",Bmass,mean,sigma1) ;

  RooRealVar  sigma2("sigma2","sigma2 of the gaussian2",0.11837);//fixed sigma
  //RooRealVar  sigma2("sigma2","sigma2 of the gaussian2",0,1);
  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 meanBkg1("meanBkg1","mean of 1st Gaussian",5.1, 5.2);
  RooRealVar sigmaBkg1("sigmaBkg1","sigma of 1st Gaussian", 0., 1.);
  RooGaussian bkgGau1("bkgGau1","bkg gaussian1 component", Bmass, meanBkg1, sigmaBkg1);

  RooRealVar nsigG("sigG","signal events", 0, 2000);
  RooRealVar nbkgExpo("nbkgExpo","number of bkg events", 0, 2000);
  RooRealVar nbkgG1("nbkgG1","number of gaussian1 bkg events", 0, 1351);

  RooAddPdf model("model","",RooArgSet(sig,bkgExpo),RooArgSet(nsigG,nbkgExpo));
  //  RooAddPdf model("model","",RooArgSet(sig1,bkgExpo,bkgGau1),RooArgSet(nsigG,nbkgExpo,nbkgG1));
  model.Print("t");

  RooFitResult* fitres = model.fitTo(*redData, Extended(true), Save(true) );

  RooPlot *xframe = Bmass.frame(Title("B_{s} invariant mass"),Bins(20) ); 
  redData->plotOn(xframe);
  model.plotOn(xframe);

  model.plotOn(xframe, Components("sig1"), LineStyle(kDashed), LineColor(kRed), FillColor(kGreen), FillStyle(3013), DrawOption("F") );
  model.plotOn(xframe, Components("bkgExpo"), LineStyle(kDotted), LineColor(kRed) );
  //model.plotOn(xframe, Components("bkgGau1"), LineStyle(kDotted), LineColor(9) );

  //  double effsigma = sqrt(sigfrac.getVal()*sigfrac.getVal()*sigma1.getVal()*sigma1.getVal()+(1-sigfrac.getVal())*(1-sigfrac.getVal())*sigma2.getVal()*sigma2.getVal());
  //  cout<<"effective sigma "<<effsigma<<endl;

  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("nbkgExpo = %.0f #pm %.0f evts", nbkgExpo.getVal(), nbkgExpo.getError()));
  //paveText->AddText(Form("nbkgG1 = %.0f #pm %.0f evts", nbkgG1.getVal(), nbkgG1.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()));
  paveText->AddText(Form("chi square = %.5f ", chiSquare.getVal())); 
  
  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();

  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("BsToPhiMuMu_data_up_ar.png");
  

}