#ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooRealVar.h" #include "RooDataSet.h" #include "RooConstVar.h" #include "RooGaussian.h" #include "RooProdPdf.h" #include "TCanvas.h" #include "TAxis.h" #include "TH1.h" #include "RooPlot.h" #include "TFile.h" #include "TStyle.h" #include "TH1F.h" using namespace RooFit ; void fit2dn() { gROOT->LoadMacro("/Users/amansangal/Desktop/Desk_top/re_reconstruct/style_file/style.C"); //--------------------BLock1-------------------------------// // declare the required variables for fitting including the ones on which cut is to be applied RooRealVar M_dz("M_dz","M_{D^{0}}(GeV/c^{2})",1.81,1.92); // for Mdz variable RooRealVar deltam("deltam","#Delta M(GeV/c^{2})",0.14,0.15); // fro deltam variable RooDataSet* ds1 = (RooDataSet*)RooDataSet::read("signal_new.dat",RooArgSet(M_dz, deltam),"Q"); ds1->Print("V"); //----------------------BLOCK 4------------------------// //--------here start constructing the models to be used for fitting -----------// //--------------------------fitter pdf for Mdz---------------------------// RooRealVar mean("mean","mean",1.863135); //--------------asym gauss1-------------------// RooRealVar sigmal1("sigmal1","sigma1",0.009544); RooRealVar sigmar1("sigmar1","#Sigma_{r1}",0.002200); RooBifurGauss asy_gauss1("asy_gauss1","asy_gauss1",M_dz,mean,sigmal1,sigmar1); //--------------asym gauss2---------------// RooRealVar sigmal2("sigmal2","sigmal",0.01739); RooRealVar sigmar2("sigmar2","sigmar",0.02275); RooBifurGauss bigauss("bigauss","bigauss",M_dz,mean,sigmal2,sigmar2); //-----------------asym gauss3-------------// RooRealVar sigmal3("sigmal3","sigmal",0.004367); RooRealVar sigmar3("sigmar3","sigmar",0.003937); RooBifurGauss asy_gauss3("asy_gauss3","asy_gauss3",M_dz,mean,sigmal3,sigmar3); //----------------CBall----------------// RooRealVar sig_cb("sig_cb","sig_cb",0.00858); RooRealVar alpha("alpha","alpha",0.896); RooRealVar n("n","n",130); RooCBShape cball("cball","cball",M_dz,mean,sig_cb,alpha,n); //--------------------fractions---------------------// RooRealVar frac_gauss1("frac_gauss1","frac_gauss1",0.219); RooRealVar frac_gauss2("frac_gauss2","frac_gauss2",0.027); RooRealVar frac_gauss3("frac_gauss3","frac_gauss3",0.57); RooRealVar frac_cb("frac_cb","frac_cb",0.37); RooFormulaVar fr2("fr2","(1-frac_gauss1)*frac_gauss2",RooArgSet(frac_gauss1,frac_gauss2)); RooFormulaVar fr3("fr3","(1-frac_gauss1)*(1-frac_gauss2)*frac_gauss3",RooArgSet(frac_gauss1,frac_gauss2,frac_gauss3)); RooFormulaVar fr4("fr4","(1-frac_gauss1)*(1-frac_gauss2)*(1-frac_gauss3)*frac_cb",RooArgSet(frac_gauss1,frac_gauss2,frac_gauss3,frac_cb)); //-------final pdf----------// RooAddPdf dg_mdz1("dg_mdz1","dg_mdz1",RooArgList(asy_gauss1,bigauss,asy_gauss3,cball),RooArgList(frac_gauss1,fr2,fr3,fr4)) ; // 3AG+ 1CB //------------------------fitter for deltam -------------------------------// //------------------(gauss+Agauss+student-t)------------------------------// RooRealVar meand("meand","meand",0.1454084); RooRealVar sigma1d("#sigma_{1}","sigma core",0.0003073); RooRealVar sigmal1d("#sigma_{l1}","sigmal1",0.0001529); RooRealVar sigmar1d("#sigma_{r1}","sigmar1",0.0000944); RooBifurGauss asygauss1d("asygauss1d","asygauss1d",deltam,meand,sigmal1d,sigmar1d); RooRealVar mean_t("mean_t","mean_t",0.1455615); RooGaussian gauss1d("gauss1d","gauss_core1",deltam,meand,sigma1d); //------------------Student-t distribution-----------------// RooRealVar r("r","r",0.760); RooRealVar sigma_t("sigma_t","sigma_t",0.0001811);//0.0001,0,1); RooGenericPdf student_t("student_t","((exp(lgamma((r+1)/2.0)-lgamma(r/2.0)))/(sqrt((22/7)*r)*sigma_t))*pow((1+(((deltam-mean_t)/sigma_t)*((deltam-mean_t)/sigma_t))/r),-(r+1)/2.0)",RooArgList(deltam,r,mean_t,sigma_t)); RooRealVar frac_gauss1d("frac_gauss1d","frac_gauss1d",0.336); RooRealVar frac_gauss2d("frac_gauss2d","frac_gauss2d",0.653); RooRealVar frac_gauss3d("frac_gauss3d","frac_gauss3",1.00); RooFormulaVar fr2d("fr2d","(1-frac_gauss1d)*frac_gauss2d",RooArgSet(frac_gauss1d,frac_gauss2d)); RooFormulaVar fr3d("fr3d","(1-frac_gauss1d)*(1-frac_gauss2d)*frac_gauss3d",RooArgSet(frac_gauss1d,frac_gauss2d,frac_gauss3d)); RooAddPdf dg_delm("dg_delm","dg_delm",RooArgList(asygauss1d,student_t,gauss1d),RooArgList(frac_gauss1d,fr2d,fr3d)); RooRealVar nsig("nsig","nsig",200000,-2000000,2000000); /////-----------Final Pdf construction --------////////// RooProdPdf mdzdelm("mdzdelm", "mdzdelm", RooArgSet(dg_mdz1,dg_delm)); RooExtendPdf mdzdel1("mdzdel1","mdzdel1",mdzdelm,nsig); RooFitResult *r1 = mdzdel1.fitTo(*ds1, Extended(kTRUE),NumCPU(6),Save()); r1->Print(); //--------------Mdz chi square-----------------------------// TCanvas *c1 = new TCanvas("c1","",550,400); c1->Divide(1,2); RooPlot *mdz_frame = M_dz.frame(); ds1->plotOn(mdz_frame);//,Binning(50));//,Binning(15)); mdzdel1.plotOn(mdz_frame, Normalization(1.0,RooAbsReal::RelativeExpected), LineColor(kBlue), LineStyle(kSolid)); mdzdel1.paramOn(mdz_frame); double chi2_ndf = mdz_frame->chiSquare();// cout<<"chi square = "<SetFillColor(10); box->SetBorderSize(1); box->SetTextAlign(12); box->SetTextSize(0.04F); box->SetFillStyle(1001); box->SetFillColor(10); TText *text = 0; Char_t buf[30]; sprintf( buf, "#chi^{2}/ndf = %f", chi2_ndf ); text = box->AddText( buf ); mdz_frame->addObject(box) ; //------------Mdz Pull Distribution---------------------// c1->cd(2)->SetPad(0.005,0.005,0.995,0.2525); RooPlot *pull1_frame = M_dz.frame(); RooHist* hpull1 = mdz_frame->pullHist(); pull1_frame->addPlotable(hpull1,"P");// p is drawing option. //----------adding reference line in pull distribution double xmin = 1.81; double xmax = 1.92; TLine *line = new TLine(xmin,0.0,xmax,0.0); TLine *line1 = new TLine(xmin,3.0,xmax,3.0); TLine *line2 = new TLine(xmin,-3.0,xmax,-3.0); line->SetLineColor(kRed); line->SetLineWidth(3); pull1_frame->GetYaxis()->SetTitle("Pull(#sigma)"); pull1_frame->GetXaxis()->CenterTitle(); pull1_frame->Draw(); line1->SetLineColor(kRed); line2->SetLineColor(kRed); line->Draw("SAME"); line1->Draw("SAME"); line2->Draw("SAME"); c1->cd(1)->SetPad(0.005,0.2525,0.995,0.995); mdz_frame->Draw(); // removed the "HIST" option c1->Draw(); //---------//--------- ------------// TCanvas *c2 = new TCanvas("c2","",550,400); c2->Divide(1,2); c2->cd(1)->SetPad(0.005,0.2525,0.995,0.995); RooPlot *deltam_frame = deltam.frame(); ds1->plotOn(deltam_frame);//,Binning(50));//,Binning(16)); mdzdel1.plotOn(deltam_frame, Normalization(1.0,RooAbsReal::RelativeExpected), LineColor(kBlue), LineStyle(kSolid)); //--------------deltam chi-square----------------------------// double chi2_ndf1 = deltam_frame->chiSquare();// cout<<"chi square1 = "<SetFillColor(10); box1->SetBorderSize(1); box1->SetTextAlign(12); box1->SetTextSize(0.04F); box1->SetFillStyle(1001); box1->SetFillColor(10); TText *text1 = 0; Char_t buf1[30]; sprintf( buf1, "#chi^{2}/ndf = %f", chi2_ndf1 ); text1 = box1->AddText( buf1 ); deltam_frame->addObject(box1) ; deltam_frame->Draw(); // removed the "HIST" option //------------deltam Pull Distribution---------------------// c2->cd(2)->SetPad(0.005,0.005,0.995,0.2525); RooPlot *pull2_frame = deltam.frame(); RooHist* hpull2 = deltam_frame->pullHist(); pull2_frame->addPlotable(hpull2,"P");// p is drawing option. //----------adding reference line in pull distribution double xmin1 = 0.14; double xmax1 = 0.15; TLine *line11 = new TLine(xmin1,0.0,xmax1,0.0); TLine *line1d = new TLine(xmin1,3.0,xmax1,3.0); TLine *line2d = new TLine(xmin1,-3.0,xmax1,-3.0); line11->SetLineColor(kRed); line11->SetLineWidth(3); pull2_frame->GetYaxis()->SetTitle("Pull(#sigma)"); pull2_frame->GetXaxis()->CenterTitle(); pull2_frame->Draw(); line1d->SetLineColor(kRed); line2d->SetLineColor(kRed); line11->Draw("SAME"); line1d->Draw("SAME"); line2d->Draw("SAME"); c2->cd(1); c2->Draw(); }