#ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooRealVar.h" #include "RooDataSet.h" #include "RooGaussian.h" #include "RooConstVar.h" #include "RooPolynomial.h" #include "RooAddPdf.h" #include "RooProdPdf.h" #include "TCanvas.h" #include "TAxis.h" #include "RooPlot.h" using namespace RooFit ; void bidi_fit_Lcnr(RooRealVar &mrho46_gRandom,RooRealVar &mrho56_gRandom) { // Create observables // RooRealVar dMass("dmfit","Lb-Lc",3.1,3.55,"GeV/c^{2}"); // RooRealVar mrho46_gRandom("mrho46_gRandom","#rho_{46}",0,3.5,"GeV/c^{2}"); // RooRealVar mrho56_gRandom("mrho56_gRandom","#rho_{56}",0,3.5,"GeV/c^{2}"); // RooRealVar m3pi("m3pi","m3#pi",0,3.5,"GeV/c^{2}"); TFile *data_file =TFile::Open("Projections_Lcnr.root"); TCanvas *canvas_2Da = new TCanvas("c2Da","2D histo Lcnr",1); TFile *data_file4 =TFile::Open("Dalitz_Rho_hepg_trig_check_mc_Lcnr.root"); TH2D *histo_data = dynamic_cast(data_file4->Get("histo1vshisto2_nocut")); // TH2* dh = histo1vshisto2_cut->createHistogram("dh",mrho46_gRandom,Binning(100),YVar(mrho56_gRandom,Binning(100))); // for data // RooDataHist dataSet_histo("dataSet_histo","dataSet_histo with mrho46_gRandom vs mrho56_gRandom",RooArgList(mrho46_gRandom,mrho56_gRandom)); RooDataHist dataSet_histo("dataSet_histo","dataSet_histo with mrho46_gRandom vs mrho56_gRandom",RooArgList(mrho46_gRandom,mrho56_gRandom)); histo_data->Draw("contz"); // histo_data->Draw("LEGO"); // dh->Draw("LEGO"); TH1 *ha = dynamic_cast(histo_data->ProjectionX("ha",1,100,"proj X")); RooDataHist dataSet_a("dataSet_a","dataSet_a",RooArgList(mrho46_gRandom),ha); TH1 *hb = dynamic_cast(histo_data->ProjectionY("hb",1,100,"proj Y")); RooDataHist dataSet_b("dataSet_b","dataSet_b",RooArgList(mrho56_gRandom),hb); TCanvas *canvas_1Da = new TCanvas("c1Da","1D histo Lcnr",1); canvas_1Da->Divide(1,2); canvas_1Da->cd(1); ha->SetTitle("ProjectionX total"); ha->Draw(); canvas_1Da->cd(2); hb->SetTitle("ProjectionY total"); hb->Draw(); TH1 *h1 = dynamic_cast(histo_data->ProjectionX("h1",1,20,"proj X")); RooDataHist dataSet1("dataSet1","dataSet1",RooArgList(mrho46_gRandom),h1); TH1 *h2 = dynamic_cast(histo_data->ProjectionY("h2",1,20,"proj Y")); RooDataHist dataSet2("dataSet2","dataSet2",RooArgList(mrho56_gRandom),h2); TH1 *h3 = dynamic_cast(histo_data->ProjectionX("h3",21,40,"proj X")); RooDataHist dataSet3("dataSet3","dataSet3",RooArgList(mrho46_gRandom),h3); TH1 *h4 = dynamic_cast(histo_data->ProjectionY("h4",21,40,"proj Y")); RooDataHist dataSet4("dataSet4","dataSet4",RooArgList(mrho56_gRandom),h4); TH1 *h5 = dynamic_cast(histo_data->ProjectionX("h5",41,60,"proj X")); RooDataHist dataSet5("dataSet5","dataSet5",RooArgList(mrho46_gRandom),h5); TH1 *h6 = dynamic_cast(histo_data->ProjectionY("h6",41,60,"proj Y")); RooDataHist dataSet6("dataSet6","dataSet6",RooArgList(mrho56_gRandom),h5); TH1 *h7 = dynamic_cast(histo_data->ProjectionX("h7",61,80,"proj X")); RooDataHist dataSet7("dataSet7","dataSet7",RooArgList(mrho46_gRandom),h7); TH1 *h8 = dynamic_cast(histo_data->ProjectionY("h8",61,80,"proj Y")); RooDataHist dataSet8("dataSet8","dataSet8",RooArgList(mrho56_gRandom),h8); TH1 *h9 = dynamic_cast(histo_data->ProjectionX("h9",81,100,"proj X")); RooDataHist dataSet9("dataSet9","dataSet9",RooArgList(mrho46_gRandom),h9); TH1 *h10 = dynamic_cast(histo_data->ProjectionY("h8",81,100,"proj Y")); RooDataHist dataSet10("dataSet10","dataSet10",RooArgList(mrho56_gRandom),h10); TCanvas *canvas_1Da = new TCanvas("c1Da","1D histo Lcnr",1); canvas_1Da->Divide(6,2); canvas_1Da->cd(1); h1->SetTitle("ProjectionX bin 1-20"); h1->Draw(); canvas_1Da->cd(2); h3->SetTitle("ProjectionX bin 21-40"); h3->Draw(); canvas_1Da->cd(3); h5->SetTitle("ProjectionX bin 41-60"); h5->Draw(); canvas_1Da->cd(4); h7->SetTitle("ProjectionX bin 61-80"); h7->Draw(); canvas_1Da->cd(5); h9->SetTitle("ProjectionX bin 81-100"); h9->Draw(); canvas_1Da->cd(6); h3->SetTitle("ProjectionX superimposed"); h7->SetLineColor(kRed); h7->DrawNormalized(); h5->SetLineColor(kBlue); h5->DrawNormalized("same"); h1->SetLineColor(kViolet); h1->DrawNormalized("same"); h3->SetLineColor(kGreen); h3->DrawNormalized("same"); h9->SetLineColor(kOrange); h9->DrawNormalized("same"); canvas_1Da->cd(7); h2->SetTitle("ProjectionY bin 1-20"); h2->Draw(); canvas_1Da->cd(8); h4->SetTitle("ProjectionY bin 21-40"); h4->Draw(); canvas_1Da->cd(9); h6->SetTitle("ProjectionY bin 41-60"); h6->Draw(); canvas_1Da->cd(10); h8->SetTitle("ProjectionY bin 61-80"); h8->Draw(); canvas_1Da->cd(11); h10->SetTitle("ProjectionY bin 81-100"); h10->Draw(); canvas_1Da->cd(12); h4->SetTitle("ProjectionY superimposed"); h6->SetLineColor(kRed); h6->DrawNormalized(); h10->SetLineColor(kBlue); h10->DrawNormalized("same"); h2->SetLineColor(kViolet); h2->DrawNormalized("same"); h8->SetLineColor(kGreen); h8->DrawNormalized("same"); h4->SetLineColor(kOrange); h4->DrawNormalized("same"); RooRealVar dm("dm","Mass shift",0,-.2,.2); RooRealVar k("k","sigma shift",1,.5,1.1); RooRealVar sig_mean7("sig_mean7","mean of the signal gaussian mass PDF",0.75,0.74,0.78); RooAddition effmean7("effmean7","mean+dm",RooArgSet(sig_mean7,dm)); RooRealVar res_sigma7("res_sigma7","sigma of the resolution function",0.23,0.22,0.26); RooProduct effsigma7("effsigma7","effsigma0",RooArgSet(res_sigma7,k)); RooGaussian gaussian7("gaussian7", "gaussian for signal peak",mrho46_gRandom,effmean7,res_sigma7); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean8("sig_mean8","mean of the signal gaussian mass PDF",1.16,1.04,1.25); RooRealVar res_sigma8("res_sigma8","sigma of the resolution function",0.63,0.6,0.7); RooProduct effsigma8("effsigma8","effsigma8",RooArgSet(res_sigma8,k)); RooAddition effmean8("effmean8","sig_mean0+dm",RooArgSet(sig_mean8,dm)); RooGaussian gaussian8("gaussian8","voigtian pdf",mrho46_gRandom,effmean8,res_sigma8); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean9("sig_mean9","mean of the signal gaussian mass PDF",1.85,1.84,1.9); RooAddition effmean9("effmean9","mean+dm",RooArgSet(sig_mean9,dm)); RooRealVar res_sigma9("res_sigma9","sigma of the resolution function",0.42,0.4,0.5); RooProduct effsigma9("effsigma9","effsigma2",RooArgSet(res_sigma9,k)); RooGaussian gaussian9("gaussian9", "gaussian for signal peak",mrho46_gRandom,effmean9,res_sigma9); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar frac2a("frac2a","frac2a",0.2,0,1); RooAddPdf model5a("model5","model5",RooArgList(gaussian7,gaussian8),frac2a); RooRealVar frac3a("frac3a","frac3a",0.2,0,1); RooAddPdf model6a("model6a","model6a",RooArgList(gaussian7,gaussian8),frac3a); dm.setConstant(); // k.setConstant(); RooFitResult *fit_res1 = model6a.fitTo(dataSet_a, Save(kTRUE), Range("SignalRegion")); TCanvas *canvas_slice_Fita = new TCanvas("canvas_sliceX_Fit1a","sliceX_Fit1",1); RooPlot *fit_frame = mrho46_gRandom.frame(Title("CDF Preliminary 2.4 fb^{-1}"),Name("fit_frame"),Bins(100)); dataSet_a.plotOn(fit_frame,Name("data")); model6a.plotOn(fit_frame,Components("model6a"),LineColor(kBlue),LineStyle(1),Name("Model")); model6a.plotOn(fit_frame,Components("gaussian7"),LineColor(kOrange),LineStyle(1),Name("gauss7")); model6a.plotOn(fit_frame,Components("gaussian8"),LineColor(kRed),LineStyle(1),Name("gauss8")); model6a.plotOn(fit_frame,Components("gaussian9"),LineColor(kGreen),LineStyle(1),Name("gauss9")); model6a.paramOn(fit_frame,Format("NELU"), Layout(0.62,1.,1.)) ; fit_frame->Draw(); sig_mean7.setConstant(); res_sigma7.setConstant(); sig_mean8.setConstant(); res_sigma8.setConstant(); sig_mean9.setConstant(); res_sigma9.setConstant(); // frac2a.setConstant(); frac3a.setConstant(); RooRealVar sig_mean7a("sig_mean7a","mean of the signal gaussian mass PDF",0.75,0.74,0.78); RooAddition effmean7a("effmean7a","mean+dm",RooArgSet(sig_mean7a,dm)); RooRealVar res_sigma7a("res_sigma7a","sigma of the resolution function",0.23,0.22,0.26); RooProduct effsigma7a("effsigma7a","effsigma7a",RooArgSet(res_sigma7a,k)); RooGaussian gaussian7a("gaussian7a", "gaussian for signal peak",mrho56_gRandom,effmean7a,res_sigma7a); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean8a("sig_mean8a","mean of the signal gaussian mass PDF",1.16,1.04,1.25); RooRealVar res_sigma8a("res_sigma8a","sigma of the resolution function",0.63,0.6,0.7); RooProduct effsigma8a("effsigma8a","effsigma8a",RooArgSet(res_sigma8a,k)); RooAddition effmean8a("effmean8a","sig_mean8a+dm",RooArgSet(sig_mean8a,dm)); RooGaussian gaussian8a("gaussian8a","voigtian pdf",mrho56_gRandom,effmean8a,res_sigma8a); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean9a("sig_mean9a","mean of the signal gaussian mass PDF",1.85,1.84,1.9); RooAddition effmean9a("effmean9a","mean+dm",RooArgSet(sig_mean9a,dm)); RooRealVar res_sigma9a("res_sigma9a","sigma of the resolution function",0.42,0.4,0.5); RooProduct effsigma9a("effsigma9a","effsigma9a",RooArgSet(res_sigma9a,k)); RooGaussian gaussian9a("gaussian9a", "gaussian for signal peak",mrho56_gRandom,effmean9a,res_sigma9a); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar frac2b("frac2b","frac2b",0.2,0,1); RooAddPdf model5b("model5b","model5b",RooArgList(gaussian7a,gaussian8a),frac2b); RooRealVar frac3b("frac3b","frac3b",0.2,0,1); RooAddPdf model6b("model6b","model6b",RooArgList(gaussian7a,gaussian8a),frac3b); dm.setConstant(); // k.setConstant(); RooFitResult *fit_res2 = model6b.fitTo(dataSet_b, Save(kTRUE), Range("SignalRegion")); TCanvas *canvas_slice_Fit2a = new TCanvas("canvas_sliceY_Fit2a","sliceY_Fit2",1); RooPlot *fit_frame2 = mrho56_gRandom.frame(Title("CDF Preliminary 2.4 fb^{-1}"),Name("fit_frame"),Bins(100)); dataSet_b.plotOn(fit_frame2,Name("data")); model6b.plotOn(fit_frame2,Components("model6b"),LineColor(kBlue),LineStyle(1),Name("Model")); model6b.plotOn(fit_frame2,Components("gaussian7a"),LineColor(kRed),LineStyle(1),Name("gauss7a")); model6b.plotOn(fit_frame2,Components("gaussian8a"),LineColor(kGreen),LineStyle(1),Name("gauss8a")); model6b.plotOn(fit_frame2,Components("gaussian9a"),LineColor(kMagenta),LineStyle(1),Name("gauss9a")); model6b.paramOn(fit_frame2,Format("NELU"), Layout(0.62,1.,1.)) ; fit_frame2->Draw(); sig_mean7a.setConstant(); res_sigma7a.setConstant(); sig_mean8a.setConstant(); res_sigma8a.setConstant(); sig_mean9a.setConstant(); res_sigma9a.setConstant(); // frac2a.setConstant(); frac3b.setConstant(); //=======================bidimensional fit======================================================== RooProdPdf totalPdf3("totalPdf3","totalPdf",RooArgSet(model6a,model6b)) ; TFile *data_file = TFile::Open("Dalitz_Rho_hepg_trig_check_mc_Lcnr.root"); TTree *data_tree = dynamic_cast(data_file->Get("lb;1")); RooDataSet data3("data3","Dataset",data_tree,RooArgSet(mrho46_gRandom,mrho56_gRandom));//,m3pi // data->reduce(Cut("m3pi>2")) ; cout << "=== Data loaded" << endl; RooFitResult *fit= totalPdf3.fitTo(data3,Extended(kFALSE),SumW2Error(kFALSE),Hesse(kTRUE),Save(kTRUE),Minos(kTRUE)); fit->Print("v"); TH2D* pdf_tot = totalPdf3.createHistogram("pdf_tot",mrho46_gRandom,YVar(mrho56_gRandom)) ; // for pdf pdf_tot->SetLineColor(kBlue) ; // TH2* dh2 = data->createHistogram("dg2",mrho46_gRandom,Binning(100),YVar(mrho56_gRandom,Binning(100))); // for data TH2* dh2 = data3.createHistogram("dh2",mrho46_gRandom,Binning(100),YVar(mrho56_gRandom,Binning(100))); // for data ,Cut("m3pi>2") TCanvas *canvas_plot2Da= new TCanvas("cplot2D","Plot",1); canvas_plot2Da->Divide(1,2); canvas_plot2Da->cd(1); pdf_tot->Draw("SURFZ"); canvas_plot2Da->cd(2); dh2->Draw("LEGO"); // dh2->Draw("same"); // canvas_plot2D->cd(3); // dh2->Draw("LEGO"); // pdf_tot->Draw("same"); // TCanvas *canvas_fitXa= new TCanvas("cfitxa","Bidimensional fit X Projection",1); RooPlot* framex = mrho46_gRandom.frame() ; data3.plotOn(framex) ; //,Cut("m3pi>2") totalPdf3.plotOn(framex) ; framex->Draw(); TCanvas *canvas_fitYa= new TCanvas("cfitya","Bidimensional fit Y Projection",1); RooPlot* framey = mrho56_gRandom.frame() ; data3.plotOn(framey) ; //,Cut("m3pi>2") totalPdf3.plotOn(framey) ; framey->Draw(); RooWorkspace *ws = new RooWorkspace("Lcnr_ws","Lcnr template WS"); ws->import(totalPdf3,RenameAllVariablesExcept("LcnrVar","mrho46_gRandom,mrho56_gRandom,dm"),RenameAllNodes("totalPdf_LcnrPdf")); ws->writeToFile("Lcnr_WS_template2D.root"); }