#ifndef __CINT__A #include "RooGlobalFunc.h" #endif #include #include #include #include #include #include #include //#include #include "RooFit.h"// #include "RooRealVar.h" #include "RooAbsPdf.h" #include "RooAddPdf.h" #include "RooDataSet.h" #include "RooDataHist.h" #include "RooGaussian.h" #include "RooGlobalFunc.h" #include "RooBifurGauss.h" #include "RooLandau.h" #include "RooChebychev.h" #include "RooFitResult.h" #include "RooNLLVar.h" #include "RooMinuit.h" #include "RooPlot.h" #include "RooArgList.h" #include "RooHist.h" #include "TLine.h" #include "TFile.h" #include "TCanvas.h" #include "TTree.h" #include "TH1.h" #include "TH1F.h" #include "TRandom.h" #include "TString.h" #include "TApplication.h" #include "TROOT.h" #include "TPad.h" using namespace RooFit ; using namespace std ; void Signal_2d() { float frac1,frac2,frac3; gROOT->LoadMacro("/Users/amansangal/Desktop/Desk_top/re_reconstruct/style_file/style.C"); // ofstream keff; // keff.open("fir_result.log"); // for (int i =1 ;i<3 ; i++) // { //--------------------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"); RooDataSet* scf = (RooDataSet*)RooDataSet::read("scf_new.dat",RooArgSet(M_dz, deltam),"Q"); ds1->append(*scf); RooDataSet* fsr = (RooDataSet*)RooDataSet::read("fsr_new.dat",RooArgSet(M_dz, deltam),"Q"); ds1->append(*fsr); //----------------------BLOCK 4------------------------// ********( !!!! very important part) //--------here start constructing the models to be used for fitting -----------// // Note: INitial value of a variable to be used in fitting is very important (Hint: there is a difference between local minimum and true(global) minimum) //-------------------------------Mdz pdf---------------------------------// //---------------- asym gauss1 -----------------------------// RooRealVar mean("mean","mean",1.86,1.85,1.87); RooRealVar sigmal1("sigmal1","sigma1",0.01,0,0.06) ; RooRealVar sigmar1("sigmar1","#Sigma_{r1}",0.02,0.001,0.06); RooBifurGauss asy_gauss1("asy_gauss1","asy_gauss1",M_dz,mean,sigmal1,sigmar1); RooRealVar frac_asymgauss2("frac_asymgauss2","frac_gauss2",0.1,0,1); RooRealVar frac_asymgauss1("frac_asymgauss1","frac_asymgauss1",0.40,0,1); //-----------------asym gauss2 -----------------------------// RooRealVar sigmal2("sigmal2","sigmal",0.015,0,0.10); RooRealVar sigmar2("sigmar2","sigmar",0.014,0,0.10); RooBifurGauss bigauss("bigauss","bigauss",M_dz,mean,sigmal2,sigmar2);// asymmetric gauss 2 //-----------------asym gauss3-------------// RooRealVar sigmal3("sigmal3","sigmal",0.0077,0,0.20);//0.0077,0,0.10 RooRealVar sigmar3("sigmar3","sigmar",0.0048,0,0.10); RooBifurGauss asy_gauss3("asy_gauss3","asy_gauss3",M_dz,mean,sigmal3,sigmar3); RooRealVar frac_asymgauss3("frac_asymgauss3","frac_asymgauss3",0.2,0,1.0); //--------------gauss1--------------// RooRealVar sigma1("sigma1","sigma1",0.001,0,0.10); RooGaussian gauss1("gauss1","gauss1",M_dz,mean,sigma1); //-----------gauss2-------------// RooRealVar sigma2("sigma2","sigma2",0.002,0,0.10); RooGaussian gauss2("gauss2","gauss1",M_dz,mean,sigma2); //-----------gauss3--------------// RooRealVar sigma3("sigma3","sigma3",0.01,0,0.10); RooGaussian gauss3("gauss3","gauss3",M_dz,mean,sigma3); //------------------Student-t distribution-----------------// RooRealVar r("r","r",2,0,10); // degree of freedom parameter for student-t distribution // RooRealVar meant("meant","meant",0,0.144,0.147); RooRealVar sigmat("sigmat","sigmat",0.01,0,0.10); //RooGenericPdf studentt("studentt","((exp(lgamma((r+1)/2.0)-lgamma(r/2.0)))/(sqrt((3.14159265359)*r)*sigmat))*pow((1+(((deltam-meant)/sigmat)*((deltam-meant)/sigmat))/r),-(r+1)/2.0)",RooArgList(deltam,r,meant,sigmat)); RooGenericPdf studentt("studentt","((exp(lgamma((r+1)/2.0)-lgamma(r/2.0)))/(sqrt((3.14159265359)*r)*sigmat))*pow((1+(((M_dz-mean)/sigmat)*((M_dz-mean)/sigmat))/r),-(r+1)/2.0)",RooArgList(M_dz,r,mean,sigmat)); //---------fractions-----------------// RooFormulaVar fr2("fr2","(1-frac_asymgauss1)*frac_asymgauss2",RooArgSet(frac_asymgauss1,frac_asymgauss2)); RooFormulaVar fr3("fr3","(1-frac_asymgauss1)*(1-frac_asymgauss2)*frac_asymgauss3",RooArgSet(frac_asymgauss1,frac_asymgauss2,frac_asymgauss3)); // RooAddPdf dg_mdz1("dg_mdz1","dg_mdz1",RooArgList(gauss1,studentt),RooArgList(frac_asymgauss1)) ; //1 Gauss + 1 student-t RooAddPdf dg_mdz1("dg_mdz1","dg_mdz1",RooArgList(asy_gauss1,bigauss,asy_gauss3),RooArgList(frac_asymgauss1,fr2)) ; //------------------------------deltam -------------------------------------// RooRealVar meand("meand","meand",0.1455,0.144,0.147); // significance of range of a variable. RooRealVar frac_gauss1d("frac_gauss1d","frac_core",0.20,0,1); RooRealVar frac_gauss2d("frac_gauss2d","frac_gauss2",0.10,0,1); // fraction for gauss 2 /////////-----------Asygauss PDF------------///////////// RooRealVar sigmald("#sigma_{ld}","sigmal",0.001,0,0.05); RooRealVar sigmard("#sigma_{rd}","sigmar",0.001,0,1); RooBifurGauss asygaussd("asygaussd","asygaussd",deltam,meand,sigmald,sigmard); RooRealVar sigmal1d("#sigma_{l1d}","sigmal1",0.0001,0,0.05); RooRealVar sigmar1d("#sigma_{r1d}","sigmar1",0.01,0,1); RooBifurGauss asygauss1d("asygauss1d","asygauss1",deltam,meand,sigmal1d,sigmar1d); //------------------Student-t distribution-----------------// RooRealVar rd("rd","rd",2.15,0,3); // degree of freedom parameter for student-t distribution RooRealVar meantd("meantd","meantd",0.1455,0.144,0.147); RooRealVar sigmatd("sigmatd","sigmatd",0.001,0,0.005); RooGenericPdf studenttd("studenttd","((exp(lgamma((rd+1)/2.0)-lgamma(rd/2.0)))/(sqrt((3.14159265359)*rd)*sigmatd))*pow((1+(((deltam-meand)/sigmatd)*((deltam-meand)/sigmatd))/rd),-(rd+1)/2.0)",RooArgList(deltam,rd,meand,sigmatd)); RooFormulaVar fr2d("fr2d","(1-frac_gauss1d)*frac_gauss2d",RooArgSet(frac_gauss1d,frac_gauss2d)); RooAddPdf dg_delm1("dg_delm1","dg_delm1",RooArgList(studenttd,asygauss1d,asygaussd),RooArgList(frac_gauss1d,fr2d)); //1student-t +2AG //----------------------final 2d pdf---------------// //RooAddPdf df_mgz("dg_mdz","dg_mdz",RooArgList(gauss1,gauss2,bigauss),RooArgList(n1,n2,n3)); unsigned long nentries = (long)ds1->numEntries(); RooRealVar nsig("nsig","nsig",500,0,1000000); //RooAddPdf dg_mdz("dz_mdz","dg_mdz",RooArgList(dg_mdz1),RooArgList(nsig)); // RooExtendPdf dg_mdz("dg_mdz","dg_mdz",dg_mdz1,nsig); // RooExtendPdf dg_delm("dg_delm","dg_delm",dg_delm1,nsig); RooProdPdf mdzdelm("mdzdelm","mdzdelm",RooArgList(dg_mdz1,dg_delm1)); RooAddPdf mdzdel1("mdzdel1","mdzdel1",mdzdelm,nsig); //----------------First call to MINUIT (no hesse and minos)----------------// mdzdel1.fitTo(*ds1,Extended(kTRUE),Minos(0),Hesse(0),NumCPU(4)); cout<<"after first minuit call"<Print("V"); //----------------First call to MINUIT ( minos for error)----------------// RooFitResult* fitres =mdzdel1.fitTo(*ds1,Extended(kTRUE),Minos(1),NumCPU(4),Save()); cout<<"after second minuit call"<Print("V"); fitres->Print("v"); //----------------getting # of parameters-----------------// //-------------getting # of degrees of freedom-----------// RooArgSet observables(M_dz,deltam); RooArgSet *flparams = mdzdel1.getParameters(observables); int nparam = (flparams->selectByAttrib("Constant",kFALSE))->getSize(); cout<<" ndf are "<Divide(1,2); RooPlot *mdz_frame = M_dz.frame(); ds1->plotOn(mdz_frame);//,Binning(50));//,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(nparam); 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); //gPad->SetLeftMargin(0.15); //pull1_frame->GetYaxis()->SetTitleOffset(1.45); pull1_frame->GetYaxis()->SetTitle("Pull(#sigma)"); pull1_frame->GetXaxis()->CenterTitle(); pull1_frame->Draw(); line1->SetLineColor(kRed); line2->SetLineColor(kRed); line1->SetLineStyle(2); line2->SetLineStyle(2); //pull1_frame->Draw(); //line->SetLineColor(kRed); line->Draw("SAME"); //pull1_frame->GetYaxis()->SetRangeUser(-5.0, 5.0); line1->Draw("SAME"); line2->Draw("SAME"); //pull1_frame->Draw(); c1->cd(1)->SetPad(0.005,0.2525,0.995,0.995); mdz_frame->Draw(); // removed the "HIST" option c1->Draw(); //-----------------------------------Deltam M frame-------------// 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_frame->Draw("HIST"); //--------------deltam chi-square----------------------------// double chi2_ndf1 = deltam_frame->chiSquare(nparam);// 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(); //------------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); //gPad->SetLeftMargin(0.15); //pull1_frame->GetYaxis()->SetTitleOffset(1.45); pull2_frame->GetYaxis()->SetTitle("Pull(#sigma)"); pull2_frame->GetXaxis()->CenterTitle(); pull2_frame->Draw(); line1d->SetLineColor(kRed); line2d->SetLineColor(kRed); line1d->SetLineStyle(2); line2d->SetLineStyle(2); //pull1_frame->Draw(); //line->SetLineColor(kRed); line11->Draw("SAME"); //pull1_frame->GetYaxis()->SetRangeUser(-5.0, 5.0); line1d->Draw("SAME"); line2d->Draw("SAME"); c2->cd(1); c2->Draw(); }