#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 part4_all1() { //--------------------BLock1-------------------------------// // declare the required variables for fitting 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 int nsig_t,ncomb_t,nrand_t,nbrcomb_t; RooDataSet* ds1 = (RooDataSet*)RooDataSet::read("signal_new.dat",RooArgSet(M_dz, deltam),"Q"); ds1->Print("V"); nsig_t = ds1->numEntries(); //----------------deltam left sidenand---------// M_dz.setRange("sbdml",1.81,1.92); deltam.setRange("sbdml",0.14,0.143); //----------------deltam right sidenand---------// M_dz.setRange("sbdmr",1.81,1.92); deltam.setRange("sbdmr",0.148,0.15); //------------mdz left side-band----------// M_dz.setRange("sbml",1.81,1.85); deltam.setRange("sbml",0.14,0.15); //------------mdz left side-band----------// M_dz.setRange("sbmr",1.88,1.92); deltam.setRange("sbmr",0.14,0.15); //-------------signal region---------// M_dz.setRange("sr",1.81,1.92); deltam.setRange("sr",0.14,0.15); //----------------------BLOCK 4------------------------// ********( !!!! very important part) //--------here start constructing the models to be used for fitting -----------// // ------1> Signal Mdz and deltam-------------// //--------------------------fitter pdf for Mdz Signal---------------------------// //-------------------( 3 Asymmetric Gauss ) -------------------// //------------------------AG1--------------------// RooRealVar mean_sig("mean_sig","mean_sig",1.86,1.85,1.87); //1.864754 RooRealVar sigmal1_sig("sigmal1_sig","sigma1_sig",0.002747);//0.01,0,0.06) ;//0.002749); RooRealVar sigmar1_sig("sigmar1_sig","#Sigma_{r1}^{sig}",0.002794);//0.02,0.001,0.06);//0.002796); RooBifurGauss asy_gauss1_sig("asy_gauss1_sig","asy_gauss1_sig",M_dz,mean_sig,sigmal1_sig,sigmar1_sig); //--------------------------AG2----------------// RooRealVar sigmal2_sig("sigmal2_sig","sigmal_sig",0.02561);//0.015,0,0.10);//0.02560); RooRealVar sigmar2_sig("sigmar2_sig","sigmar_sig",0.02033);//0.014,0,0.10);//0.02032); RooBifurGauss bigauss_sig("bigauss_sig","bigauss_sig",M_dz,mean_sig,sigmal2_sig,sigmar2_sig);// asymmetric gauss 2 //-----------------asym gauss3-------------// RooRealVar sigmal3_sig("sigmal3_sig","sigmal_sig",0.00606);//0.0077,0,0.20);//0.00607); RooRealVar sigmar3_sig("sigmar3_sig","sigmar_sig",0.00617);//0.0048,0,0.10); RooBifurGauss asy_gauss3_sig("asy_gauss3_sig","asy_gauss3_sig",M_dz,mean_sig,sigmal3_sig,sigmar3_sig); //---------fractions-----------------// RooRealVar frac_asymgauss2_sig("frac_asymgauss2_sig","frac_gauss2_sig",0.1312);//0.1,0,0.65); RooRealVar frac_asymgauss1_sig("frac_asymgauss1_sig","frac_asymgauss1_sig",0.679);//0.35,0,1); RooFormulaVar fr2_sig("fr2_sig","(1-frac_asymgauss1_sig)*frac_asymgauss2_sig",RooArgList(frac_asymgauss1_sig,frac_asymgauss2_sig)); //--------------final pdf for Mdz signal------------------// RooAddPdf dg_mdz1_sig("dg_mdz1_sig","dg_mdz1_sig",RooArgList(asy_gauss1_sig,bigauss_sig,asy_gauss3_sig),RooArgList(frac_asymgauss1_sig,fr2_sig)) ; // 3 Asy Gauss cout<<"one"< random pi slow background Mdz and deltam-------------// //------------------------fitter for deltam random pi slow -------------------------------// //-------------------------Threshold Function----------------------// RooRealVar pip_mass_rand("pip_mass_rand","#pi^{+} Mass",0.13957061); // pdg value for pip mass is 0.13957018 RooRealVar alpha1_rand("alpha1_rand","alpha1_rand",3.6);//-11,-35,100); RooGenericPdf dg_delm_rand("dg_delm_rand","dg_delm_rand","pow(deltam-0.13957061,0.5)+alpha1_rand*pow((deltam-0.13957061),1.5)",RooArgSet(deltam,alpha1_rand)); // ------3> combinatoric background Mdz and deltam-------------// //--------------------------fitter pdf for Mdz , combinatorial background---------------------------// //-------------------Chebychev polynomial----------------------// RooRealVar a1_comb("a1_comb","a1_comb",1,0,5);//0.173); RooRealVar a2_comb("a2_comb","a2_comb",0.5,-1,1);//-0.0309); //RooRealVar a3("a3","a3",0.1,-10,10); RooChebychev dg_mdz1_comb("dg_mdz1_comb","dg_mdz1_comb",M_dz,RooArgList(a1_comb,a2_comb)); //------------------------fitter for deltam, combinatorial background -------------------------------// //-------------------Threshold Function----------------------// RooRealVar pip_mass_comb("pip_mass_comb","#pi^{+} Mass",0.13957061); // pdg value for pip mass is 0.13957018 RooRealVar alpha1_comb("alpha1_comb","alpha1_comb",-38.41);//-20,-50,20);//-38.41); RooGenericPdf threshold_comb("threshold_comb","threshold_comb","pow(deltam-0.13957061,0.5)+alpha1_comb*pow((deltam-0.13957061),1.5)",RooArgSet(deltam,alpha1_comb)); RooRealVar meand_comb("meand_comb","meand_comb",0.1455,0.14,0.15); RooRealVar sigmad_comb("sigmad_comb","sigmad_comb",0.01,0,0.10); RooGaussian gauss_comb("gauss_comb","gauss_comb",deltam,meand_comb,sigmad_comb); RooRealVar fracg_comb("fracg_comb","",0.2,0,1); RooAddPdf dg_delm_comb("dg_delm_comb","dg_delm_comb",RooArgList(threshold_comb,gauss_comb),RooArgList(fracg_comb)); //-------------------Gaussian Function ----------------------// /* RooRealVar delm1_comb("delm1_comb","delm1_comb",0.14,0.15); // to change the range of gaussian RooRealVar mean_comb("mean_comb","mean_comb",0.14506); RooRealVar sigma_comb("#sigma_comb","#sigma_comb",0.00100); RooGaussian gauss1d_comb("gauss1d_comb","gauss1d_comb",deltam,mean_comb,sigma_comb); RooRealVar frac_threshold_comb("frac_{threshold}^{comb}","frac_{threshold}",0.9809); //---------------------Adding the two: Threshold + Gaussian----------------------------------// RooAddPdf dg_delm_comb("dg_delm_comb","dg_delm_comb",RooArgList(threshold_comb,gauss1d_comb),RooArgList(frac_threshold_comb));*/ //------------------3ks bkg-------------- // //--------------------------fitter pdf for Mdz---------------------------// //-----------------------Gauss pdf -------------------------// //---------------- gauss1 -----------------------------// RooRealVar mean_3ks("mean_3ks","mean_3ks",1.85711);//1.86,1.85,1.87); RooRealVar sigma1_3ks("sigma1_3ks","sigma1_3ks",0.02063);//0.01,0,0.10); RooGaussian gauss1_3ks("gauss1_3ks","gauss1_3ks",M_dz,mean_3ks,sigma1_3ks); //------------------------fitter for deltam -------------------------------// //-------------student-t-------------// RooRealVar meand_3ks("meand_3ks","meand_3ks",0.145475);//0.1455,0.144,0.147); // significance of range of a variable. //------------------Student-t distribution-----------------// RooRealVar rd_3ks("rd_3ks","rd_3ks",1.67);//2.15,0,3); // degree of freedom parameter for student-t distribution RooRealVar sigmatd_3ks("sigmatd_3ks","sigmatd_3ks",0.000523);//0.001,0,0.005); RooGenericPdf studenttd_3ks("studenttd_3ks","((exp(lgamma((rd_3ks+1)/2.0)-lgamma(rd_3ks/2.0)))/(sqrt((3.14159265359)*rd_3ks)*sigmatd_3ks))*pow((1+(((deltam-meand_3ks)/sigmatd_3ks)*((deltam-meand_3ks)/sigmatd_3ks))/rd_3ks),-(rd_3ks+1)/2.0)",RooArgList(deltam,rd_3ks,meand_3ks,sigmatd_3ks)); //-----------------broken charm bkg------------------------// //-------mdz--------------------// RooRealVar mean_bc("mean_bc","mean_bc",1.86493);//1.86,1.85,1.87); //--------------gauss1--------------// RooRealVar sigma1_bc("sigma1_bc","sigma1_bc",0.052);//0.005,0,0.10); RooGaussian gauss1_bc("gauss1_bc","gauss1_bc",M_dz,mean_bc,sigma1_bc); //-----------gauss2-------------// RooRealVar sigma2_bc("sigma2_bc","sigma2_bc",0.00836);//0.002,0,0.10); RooGaussian gauss2_bc("gauss2_bc","gauss1_bc",M_dz,mean_bc,sigma2_bc); RooRealVar frac_asymgauss1_bc("frac_asymgauss1_bc","frac_asymgauss1_bc",0.392);//0.40,0,1); RooAddPdf dg_mdz_bc("dg_mdz_bc","dg_mdz_bc",RooArgList(gauss1_bc,gauss2_bc),RooArgList(frac_asymgauss1_bc)) ; //----------------------deltam pdf-------------// RooRealVar meand_bc("meand_bc","meand_bc",0.145422);//0.1455,0.144,0.147); // significance of range of a variable. //------------------Student-t distribution-----------------// RooRealVar rd_bc("rd_bc","rd_bc",0.600);//2.15,0,3); // degree of freedom parameter for student-t distribution RooRealVar sigmatd_bc("sigmatd_bc","sigmatd_bc",0.000216);//0.001,0,0.005); RooGenericPdf studenttd_bc("studenttd_bc","((exp(lgamma((rd_bc+1)/2.0)-lgamma(rd_bc/2.0)))/(sqrt((3.14159265359)*rd_bc)*sigmatd_bc))*pow((1+(((deltam-meand_bc)/sigmatd_bc)*((deltam-meand_bc)/sigmatd_bc))/rd_bc),-(rd_bc+1)/2.0)",RooArgList(deltam,rd_bc,meand_bc,sigmatd_bc)); //---------------------done-------------------------------------------------// RooProdPdf mdzdelm_sig("mdzdelm_sig", "mdzdelm", RooArgSet(dg_mdz1_sig,dg_delm_sig)); // 2D pdf (deltam * mdz) for signal RooProdPdf mdzdelm_rand("mdzdelm_rand", "mdzdelm", RooArgSet(dg_mdz1_sig,dg_delm_rand)); // 2D pdf (deltam * mdz) for rando //RooProdPdf mdzdelm_comb("mdzdelm_comb", "mdzdelm", RooArgSet(dg_mdz1_comb,threshold_comb)); // 2D pdf (deltam * mdz) for combina RooProdPdf mdzdelm_comb("mdzdelm_comb", "mdzdelm", RooArgSet(dg_mdz1_comb,dg_delm_comb)); RooProdPdf mdzdelm_d3k("mdzdelm_d3k", "mdzdelm", RooArgSet(gauss1_3ks,studenttd_3ks)); RooProdPdf mdzdelm_bc("mdzdelm_bc", "mdzdelm", RooArgSet(dg_mdz_bc,studenttd_bc)); RooRealVar num_sig("num_sig","num_sig",1,0,300000); RooRealVar num_comb("num_comb","num_comb", 14000,0,300000); RooRealVar num_rand("num_rand","num_rand", 0);//10, 0,1000000); RooRealVar num_d3k("num_d3k","num_d3k",0);//130);// fixed RooRealVar num_bc("num_bc","num_bc",0);//161); /* RooArgList shapes; RooArgList yields; shapes.add(mdzdelm_sig); shapes.add(mdzdelm_comb); shapes.add(mdzdelm_rand); // shapes.add(mdzdelm_d3k); yields.add(num_sig); yields.add(num_comb); yields.add(num_rand); // yields.add(num_d3k); RooAddPdf mdzdel1("mdzdel1", "", shapes, yields);*/ RooExtendPdf mdzdelm_sig_e("mdzdelm_sig_e","mdzdelm_sig_e",mdzdelm_sig,num_sig) ; RooExtendPdf mdzdelm_comb_e("mdzdelm_comb_e","mdzdelm_comb_e",mdzdelm_comb,num_comb) ; //RooExtendPdf mdzdel1("mdzdel1","mdzdel1",mdzdelm_comb,num_comb) ; // for debugging the combinatorial fit // RooAddPdf mdzdel1("mdzdel1","mdzdel1",mdzdelm_comb,num_comb) ; // for debugging the combinatorial fit RooExtendPdf mdzdelm_rand_e("mdzdelm_rand_e","mdzdelm_rand_e",mdzdelm_rand,num_rand) ; RooExtendPdf mdzdelm_d3k_e("mdzdelm_d3k_e","mdzdelm_d3k_e",mdzdelm_d3k,num_d3k) ; RooExtendPdf mdzdelm_bc_e("mdzdelm_bc_e","mdzdelm_bc_e",mdzdelm_bc,num_bc) ; /* RooExtendPdf mdzdelm_sig_e("mdzdelm_sig_e","mdzdelm_sig_e",dg_mdz1_sig,num_sig) ; RooExtendPdf mdzdelm_comb_e("mdzdelm_comb_e","mdzdelm_comb_e",dg_mdz1_comb,num_comb) ; RooExtendPdf mdzdelm_rand_e("mdzdelm_rand_e","mdzdelm_rand_e",dg_mdz1_sig,num_rand) ; RooExtendPdf mdzdelm_d3k_e("mdzdelm_d3k_e","mdzdelm_d3k_e",gauss1_3ks,num_d3k) ; RooExtendPdf mdzdelm_bc_e("mdzdelm_bc_e","mdzdelm_bc_e",dg_mdz_bc,num_bc) ;*/ // RooAddPdf mdzdel1("mdzdel1", "mdzdel1", RooArgList(mdzdelm_sig_e,mdzdelm_rand_e,mdzdelm_comb_e,mdzdelm_d3k_e,mdzdelm_bc_e)); //RooAddPdf mdzdel1("mdzdel1", "mdzdel1", RooArgList(mdzdelm_comb_e)); cout<<"before fit to"<Print("v"); // printing the fit results Int_t st = rs1->status(); cout<<"status is "<< st<covQual(); cout<<"cov quality is "<< st<selectByAttrib("Constant",kFALSE))->getSize(); cout<<" ndf are "<sumEntries("", "sbmr,sbml"); cout<<" nentries in data are "<Divide(1,2); //c1_1->SetPad(0.005,0.2525,0.995,0.995); //c1_2->SetPad(0.005,0.005,0.995,0.2525); //c1->cd(1); RooPlot *mdz_frame = M_dz.frame(); ds1->plotOn(mdz_frame);//,Binning(50));//,Binning(15));Range(“sb_lo,sb_hi”),NormRange(“sb_lo,sb_hi”) mdzdel1.plotOn(mdz_frame, Range("sbml,sbmr"),NormRange("sbml,sbmr"),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); /* mdzdel1.plotOn(mdz_frame,Name("Signal"), Components(mdzdelm_sig), LineColor(kGreen), LineStyle(8)); mdzdel1.plotOn(mdz_frame,Name("Combinatorial"), Components(mdzdelm_comb), LineColor(kRed), LineStyle(3)); mdzdel1.plotOn(mdz_frame,Name("Random #pi_{slow}"), Components(mdzdelm_rand), LineColor(kMagenta), LineStyle(5)); mdzdel1.plotOn(mdz_frame,Name("D^{0}->3K_{s}^{0}"), Components(mdzdelm_d3k), LineColor(kBlack), LineStyle(5)); mdzdel1.plotOn(mdz_frame,Name("Broken charm"), Components(mdzdelm_bc), LineColor(28), LineStyle(5));*/ mdz_frame->Draw(); // removed the "HIST" option TLegend *leg = new TLegend(0.5,0.5,0.6,0.6); leg->AddEntry("Signal","Signal","L"); leg->AddEntry("Combinatorial","Combinatorial","L"); leg->AddEntry("Random #pi_{slow}","Random #pi_{slow}","L"); leg->AddEntry("D^{0}->3K_{s}^{0}","D^{0}->3K_{s}^{0}","L"); leg->AddEntry("Broken charm","Broken charm","L"); leg->Draw(); 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, Range("sbml,sbmr"),NormRange("sbml,sbmr"), LineColor(kBlue), LineStyle(kSolid),Normalization(nData, RooAbsReal::NumEvent)); //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(); // 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); //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); /* mdzdel1.plotOn(deltam_frame,Name("Signal"), Components(mdzdelm_sig), LineColor(kGreen), LineStyle(8)); mdzdel1.plotOn(deltam_frame,Name("Combinatorial"), Components(mdzdelm_comb), LineColor(kRed), LineStyle(3)); mdzdel1.plotOn(deltam_frame,Name("Random #pi_{slow}"), Components(mdzdelm_rand), LineColor(kMagenta), LineStyle(5)); mdzdel1.plotOn(deltam_frame,Name("D^{0}->3K_{s}^{0}"), Components(mdzdelm_d3k), LineColor(kBlack), LineStyle(5)); mdzdel1.plotOn(deltam_frame,Name("Broken charm"), Components(mdzdelm_bc), LineColor(28), LineStyle(5));*/ deltam_frame->Draw(); // removed the "HIST" option TLegend *leg1 = new TLegend(0.5,0.5,0.6,0.6); leg1->AddEntry("Signal","Signal","L"); leg1->AddEntry("Combinatorial","Combinatorial","L"); leg1->AddEntry("Random #pi_{slow}","Random #pi_{slow}","L"); leg1->AddEntry("D^{0}->3K_{s}^{0}","D^{0}->3K_{s}^{0}","L"); leg1->AddEntry("Broken charm","Broken charm","L"); leg1->Draw(); c2->Draw(); //--------printing the values of residual and pull for yields----------// /* RooRealVar* nsig_fitresult = (RooRealVar*) rs1->floatParsFinal().find(num_sig); cout<<" signal yield is "<getVal()<<"+-"<getError()<getVal()-nsig_t)<<" "<<(nsig_fitresult->getVal()-nsig_t)/nsig_fitresult->getError()<floatParsFinal().find(num_comb); cout<<" combi yield is "<getVal()<<"+-"<getError()<getVal()-(ncomb_t-nbrcomb_t))<<" "<<(ncomb_fitresult->getVal()-(ncomb_t-nbrcomb_t))/ncomb_fitresult->getError()<floatParsFinal().find(num_rand); cout<<" rand yield is "<getVal()<<"+-"<getError()<getVal()-nrand_t)<<" "<<(nrand_fitresult->getVal()-nrand_t)/nrand_fitresult->getError()<