// Relativistic Breit-Weigner pdf to fit m(KK) distribution const char* file = "/eos/user/r/rraturi/f2pmumu/red_ntuples"; const char* plot = "/eos/user/r/rraturi/www/analysis/2016/tmva"; using namespace RooFit; void plot_pull(RooPlot* pull, double xmin, double xmax) { pull->SetMinimum(-5); pull->SetMaximum(5); pull->SetFillColor(kBlue); pull->GetYaxis()->SetNdivisions(5); pull->SetTitle(""); pull->GetYaxis()->SetTitle("Pull"); pull->GetYaxis()->SetTitleSize(0.9); pull->GetYaxis()->SetLabelSize(0.1); pull->GetXaxis()->SetLabelSize(0.1); pull->GetXaxis()->SetTitle("#bf{m(K^{+}K^{-}) [GeV]}"); pull->GetXaxis()->SetLimits(xmin, xmax); } void plot_lines(double xmin, double xmax) { TLine* l1 = new TLine(xmin, +3, xmax, +3); TLine* l2 = new TLine(xmin, -3, xmax, -3); l1->SetLineColor(2); l2->SetLineColor(2); TLine* l3 = new TLine(xmin, 0, xmax, 0); l3->SetLineColor(2); l1->Draw(); l2->Draw(); l3->Draw(); } void rel_bw() { gROOT->SetBatch(1); TChain *ch = new TChain("tree"); ch->Add(Form("%s/JPsif2p_MC_ForFit.root", file)); TTree *tr = ch; int nentries_ = tr->GetEntries(); cout << "\n=> total entries in signal tree = " << nentries_ << endl; double ph_min(1.353535), ph_max(1.78); double bm_min(5.27), bm_max(5.45); RooRealVar Phimass("Phimass","#bf{m(K^{+}K^{-}) [GeV]}", ph_min, ph_max); RooRealVar Bmass("Bmass", "", bm_min, bm_max); RooArgSet observables(Phimass, Bmass); RooDataSet data("data","dataset with Bmass", ch, observables); TCut c1 = Form("Phimass>%f && Phimass<%f",ph_min,ph_max); TCut cut2 = Form("Bmass>%f && Bmass<%f",bm_min,bm_max); TCut cutTotal = c1 && cut2 ; RooDataSet *redData = (RooDataSet*)data.reduce(cutTotal); RooRealVar mn("mn","common means for Crystal Balls", 1.525, ph_min, ph_max); RooRealVar sigma1("sigma1","sigma of CB1", 0.05, 0., 0.1); RooRealVar sigma2("sigma2","sigma of CB2", 0.07, 0., 0.1); RooFormulaVar gamma1("gamma1", "sqrt(mn*mn*(mn*mn+sigma1*sigma1))", RooArgList(mn, sigma1)); RooFormulaVar gamma2("gamma2", "sqrt(mn*mn*(mn*mn+sigma2*sigma2))", RooArgList(mn, sigma2)); RooFormulaVar k1("k1", "(2.83*mn*sigma1*gamma1)/(3.14*sqrt(mn*mn+gamma1))", RooArgList(mn, sigma1, gamma1)); RooFormulaVar k2("k2", "(2.83*mn*sigma2*gamma2)/(3.14*sqrt(mn*mn+gamma2))", RooArgList(mn, sigma2, gamma2)); RooGenericPdf rel_bw1("rel_bw1", "", "k1/((Phimass**2-mn**2)**2+(mn*sigma1)**2)", RooArgList(Phimass, mn, sigma1, k1)); RooGenericPdf rel_bw2("rel_bw2", "", "k2/((Phimass**2-mn**2)**2+(mn*sigma2)**2)", RooArgList(Phimass, mn, sigma2, k2)); RooRealVar frac("frac","fraction of CB", 0.4, 0.01, 1.); RooAddPdf rel_bw("rel_bw","", RooArgList(rel_bw1,rel_bw2), RooArgList(frac)); RooRealVar nsig("nsig","nsig",1E3,500,1.1*tr->GetEntries()); RooExtendPdf RBW("RBW", "RBW", rel_bw, nsig); RBW.Print("t"); TCanvas *pdf = new TCanvas("pdf", "", 800, 700); gPad->SetTickx(); gPad->SetTicky(); RooPlot *xpdf = Phimass.frame(Title(""), Bins(200)); xpdf->SetTitle(""); xpdf->GetYaxis()->SetTitleOffset(1.55); RBW.plotOn(xpdf); xpdf->Draw(); pdf->SaveAs(Form("%s/rel_bw_distribution.pdf", plot)); pdf->SaveAs(Form("%s/rel_bw_distribution.png", plot)); RooFitResult* fitres = RBW.fitTo(*redData, Extended(true), Minos(true), Save(true)); TCanvas *c = new TCanvas("c","c", 800, 700); TPad *p1 = new TPad("p1","p1", 0.01, 0.25, 0.995, 0.97); TPad *p2 = new TPad("p2","p2", 0.01, 0.02, 0.995, 0.24); p1->Draw(); p2->Draw(); p1->cd(); gPad->SetTickx(); gPad->SetTicky(); gPad->SetBottomMargin(0); RooPlot *xframe = Phimass.frame(Title(""), Bins(200)); xframe->SetTitle(""); //xframe->GetXaxis()->SetTitle("#bf{m(K^{+}K^{-}) [GeV]}"); xframe->GetYaxis()->SetTitleOffset(1.55); redData->plotOn(xframe,RooFit::Name("data")); RBW.plotOn(xframe, RooFit::Name("Full PDF"), LineColor(kBlue)); RooHist* hpull = xframe->pullHist(); printf("**************************************"); printf("\n #sigma_{1}: %.7f", sigma1.getVal()); printf("\n #sigma_{2}: %.7f", sigma2.getVal()); double eff_sigma = sqrt(frac.getVal()*pow(sigma1.getVal(),2) + (1 - frac.getVal()) * pow(sigma2.getVal(),2)); printf("\n #sigma_{eff}: %.7f", eff_sigma); printf("\n************************************\n"); int nFloatParam = fitres->floatParsFinal().getSize(); cout << "number of floating parameters-> " << nFloatParam << endl; double chi2dof = xframe->chiSquare(nFloatParam); std::cout<<"\n"<SetBorderSize(0.0); pText->SetFillColor(kWhite); pText->SetFillStyle(0); pText->SetTextSize(0.02); pText->AddText(Form("N_{sig} = %.0f #pm %.0f", nsig.getVal(), nsig.getError())); pText->AddText(Form("#mu = %.7f #pm %.7f GeV" , mn.getVal() , mn.getError())); pText->AddText(Form("#sigma_{1} = %.7f #pm %.7f GeV", sigma1.getVal(), sigma1.getError())); pText->AddText(Form("#sigma_{2} = %.7f #pm %.7f GeV", sigma2.getVal(), sigma2.getError())); pText->AddText(Form("#sigma_{eff} = %.7f GeV", eff_sigma)); pText->AddText(Form("#chi^{2}/dof = %.5f ", chi2dof )); pText->Draw(); TLatex *mark = new TLatex(); mark->SetNDC(true); double startY = 0.92; mark->SetTextFont(42); mark->SetTextSize(0.035); mark->DrawLatex(0.12,startY,"#bf{CMS} #it{Simulation}"); mark->DrawLatex(0.70,startY,"#scale[0.8]{65.9 fb^{-1} (13 TeV)}"); mark->Draw(); p2->cd(); gPad->SetTickx(); gPad->SetTicky(); gPad->SetTopMargin(0); RooPlot* pull = Phimass.frame(Title("")); pull->addPlotable(hpull, "BX"); hpull->SetFillColor(kBlue); plot_pull(pull, ph_min, ph_max); pull->GetXaxis()->SetTitle("#bf{m(K^{+}K^{-}) [GeV]}"); pull->Draw("P"); plot_lines(ph_min, ph_max); cout << "error #gamma1" << gamma1.defaultErrorLevel() << endl; cout << "error #gamma2" << gamma2.defaultErrorLevel() << endl; c->SaveAs(Form("%s/Dikaon_relBW_bdtr.pdf", plot)); c->SaveAs(Form("%s/Dikaon_relBW_bdtr.png", plot)); }