#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_Lcrho(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_Lcrho.root"); TCanvas *canvas_2Db = new TCanvas("c2Db","2D histo Lcrho",1); TFile *data_file4 =TFile::Open("Dalitz_Rho_hepg_trig_check_mc_Lcrho.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_1Db = new TCanvas("c1Db","1D histo Lcrho",1); canvas_1Db->Divide(1,2); canvas_1Db->cd(1); ha->SetTitle("ProjectionX total"); ha->Draw(); canvas_1Db->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_1Db = new TCanvas("c1Db","1D histo Lcrho",1); canvas_1Db->Divide(6,2); canvas_1Db->cd(1); h1->SetTitle("ProjectionX bin 1-20"); h1->Draw(); canvas_1Db->cd(2); h3->SetTitle("ProjectionX bin 21-40"); h3->Draw(); canvas_1Db->cd(3); h5->SetTitle("ProjectionX bin 41-60"); h5->Draw(); canvas_1Db->cd(4); h7->SetTitle("ProjectionX bin 61-80"); h7->Draw(); canvas_1Db->cd(5); h9->SetTitle("ProjectionX bin 81-100"); h9->Draw(); canvas_1Db->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_1Db->cd(7); h2->SetTitle("ProjectionY bin 1-20"); h2->Draw(); canvas_1Db->cd(8); h4->SetTitle("ProjectionY bin 21-40"); h4->Draw(); canvas_1Db->cd(9); h6->SetTitle("ProjectionY bin 41-60"); h6->Draw(); canvas_1Db->cd(10); h8->SetTitle("ProjectionY bin 61-80"); h8->Draw(); canvas_1Db->cd(11); h10->SetTitle("ProjectionY bin 81-100"); h10->Draw(); canvas_1Db->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_mean0("sig_mean0","mean of the signal gaussian mass PDF",0.52,0.5,0.55); RooAddition effmean0("effmean0","mean+dm",RooArgSet(sig_mean0,dm)); RooRealVar res_sigma0("res_sigma0","sigma of the resolution function",0.1,0.098,0.15); RooProduct effsigma0("effsigma0","effsigma0",RooArgSet(res_sigma0,k)); RooGaussian gaussian0("gaussian0", "gaussian for signal peak",mrho46_gRandom,effmean0,effsigma0); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean1("sig_mean1","mean of the signal gaussian mass PDF",0.77,0.7,0.9); RooRealVar res_sigma1("res_sigma1","sigma of the resolution function",0.066,0.05,0.075); RooProduct effsigma1("effsigma1","effsigma1",RooArgSet(res_sigma1,k)); RooAddition effmean1("effmean1","sig_mean0+dm",RooArgSet(sig_mean1,dm)); RooGaussian gaussian1("gaussian1","voigtian pdf",mrho46_gRandom,effmean1,effsigma1); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean1a("sig_mean1a","mean of the signal gaussian mass PDF",0.88,0.8,0.95); RooAddition effmean1a("effmean1a","mean+dm",RooArgSet(sig_mean1a,dm)); RooRealVar res_sigma1a("res_sigma1a","sigma of the resolution function",0.22,0.15,0.3); RooProduct effsigma1a("effsigma1a","effsigma1a",RooArgSet(res_sigma1a,k)); RooGaussian gaussian1a("gaussian1a", "gaussian for signal peak",mrho46_gRandom,effmean1a,effsigma1a); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean1b("sig_mean1b","mean of the signal gaussian mass PDF",1.37,1.3,1.45); RooAddition effmean1b("effmean1b","mean+dm",RooArgSet(sig_mean1b,dm)); RooRealVar res_sigma1b("res_sigma1b","sigma of the resolution function",0.55,0.45,0.65); RooProduct effsigma1b("effsigma1b","effsigma1",RooArgSet(res_sigma1b,k)); RooGaussian gaussian1b("gaussian1b", "gaussian for signal peak",mrho46_gRandom,effmean1b,effsigma1b); RooRealVar frac2("frac2","frac2",0.2,0,1); RooAddPdf model5("model5","model5",RooArgList(gaussian0,gaussian1),frac2); RooRealVar frac3("frac3","frac3",0.2,0,1); RooAddPdf model6("model6","model6",RooArgList(model5,gaussian1a),frac3); RooRealVar peak_frac1("peak_frac1","peak_frac1",0.2,0,1); RooAddPdf model3("model3","model3",RooArgList(model6,gaussian1b),peak_frac1); dm.setConstant(); k.setConstant(); RooFitResult *fit_res1 = model3.fitTo(dataSet_a, Save(kTRUE), Range("SignalRegion")); TCanvas *canvas_slice_Fitb = new TCanvas("canvas_sliceX_Fit1b","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")); model3.plotOn(fit_frame,Components("model3"),LineColor(kBlue),LineStyle(1),Name("Model")); model3.plotOn(fit_frame,Components("gaussian0"),LineColor(kOrange),LineStyle(1),Name("gauss0")); model3.plotOn(fit_frame,Components("gaussian1"),LineColor(kRed),LineStyle(1),Name("gauss1")); model3.plotOn(fit_frame,Components("gaussian1a"),LineColor(kGreen),LineStyle(1),Name("gauss1a")); model3.plotOn(fit_frame,Components("gaussian1b"),LineColor(kMagenta),LineStyle(1),Name("gauss1b")); model3.paramOn(fit_frame,Format("NELU"), Layout(0.62,1.,1.)) ; fit_frame->Draw(); sig_mean0.setConstant(); res_sigma0.setConstant(); sig_mean1a.setConstant(); res_sigma1a.setConstant(); sig_mean1b.setConstant(); res_sigma1b.setConstant(); sig_mean1.setConstant(); res_sigma1.setConstant(); frac2.setConstant(); frac3.setConstant(); peak_frac1.setConstant(); RooRealVar sig_mean0c("sig_mean0c","mean of the signal gaussian mass PDF",0.52,0.5,0.55); RooAddition effmean0c("effmean0c","mean+dm",RooArgSet(sig_mean0c,dm)); RooRealVar res_sigma0c("res_sigma0c","sigma of the resolution function",0.1,0.098,0.15); RooProduct effsigma0c("effsigma0c","effsigma0c",RooArgSet(res_sigma0c,k)); RooGaussian gaussian0c("gaussian0c", "gaussian for signal peak",mrho56_gRandom,effmean0c,effsigma0c); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean1c("sig_mean1c","mean of the signal gaussian mass PDF",0.77,0.7,0.9); RooRealVar res_sigma1c("res_sigma1c","sigma of the resolution function",0.066,0.05,0.075); RooProduct effsigma1c("effsigma1c","effsigma1",RooArgSet(res_sigma1c,k)); RooAddition effmean1c("effmean1","sig_mean0+dm",RooArgSet(sig_mean1c,dm)); RooGaussian gaussian1c("gaussian1c","voigtian pdf",mrho56_gRandom,effmean1c,effsigma1c); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean2c("sig_mean2c","mean of the signal gaussian mass PDF",0.88,0.8,0.95); RooAddition effmean2c("effmean2c","mean+dm",RooArgSet(sig_mean2c,dm)); RooRealVar res_sigma2c("res_sigma2c","sigma of the resolution function",0.22,0.15,0.3); RooProduct effsigma2c("effsigma2c","effsigma2",RooArgSet(res_sigma2c,k)); RooGaussian gaussian2c("gaussian2c", "gaussian for signal peak",mrho56_gRandom,effmean2c,effsigma2c); //////////////////////////////////////////////////////////////////////////////////////////////// RooRealVar sig_mean3c("sig_mean3c","mean of the signal gaussian mass PDF",1.37,1.3,1.45); RooAddition effmean3c("effmean3","mean+dm",RooArgSet(sig_mean3c,dm)); RooRealVar res_sigma3c("res_sigma3c","sigma of the resolution function",0.55,0.45,0.65); RooProduct effsigma3c("effsigma3c","effsigma3c",RooArgSet(res_sigma3c,k)); RooGaussian gaussian3c("gaussian3c", "gaussian for signal peak",mrho56_gRandom,effmean3c,effsigma3c); RooRealVar frac2c("frac2c","frac2c",0.2,0,1); RooAddPdf model5c("model5c","model5c",RooArgList(gaussian0c,gaussian1c),frac2c); RooRealVar frac3c("frac3c","frac3c",0.2,0,1); RooAddPdf model6c("model6c","model6c",RooArgList(model5c,gaussian2c),frac3c); RooRealVar peak_frac1c("peak_frac1c","peak_frac1",0.2,0,1); RooAddPdf model3c("model3c","model3c",RooArgList(model6c,gaussian3c),peak_frac1c); dm.setConstant(); k.setConstant(); RooFitResult *fit_res2 = model3c.fitTo(dataSet_b, Save(kTRUE), Range("SignalRegion")); TCanvas *canvas_slice_Fit2b= new TCanvas("canvas_sliceY_Fit2b","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")); model3c.plotOn(fit_frame2,Components("model3c"),LineColor(kBlue),LineStyle(1),Name("Model")); model3c.plotOn(fit_frame2,Components("gaussian0c"),LineColor(kRed),LineStyle(1),Name("gauss2")); model3c.plotOn(fit_frame2,Components("gaussian1c"),LineColor(kGreen),LineStyle(1),Name("gauss2a")); model3c.plotOn(fit_frame2,Components("gaussian2c"),LineColor(kMagenta),LineStyle(1),Name("gauss2a")); model3c.plotOn(fit_frame2,Components("gaussian3c"),LineColor(kOrange),LineStyle(1),Name("gauss2a")); model3c.paramOn(fit_frame2,Format("NELU"), Layout(0.62,1.,1.)) ; fit_frame2->Draw(); sig_mean0c.setConstant(); res_sigma0c.setConstant(); sig_mean1c.setConstant(); res_sigma1c.setConstant(); sig_mean2c.setConstant(); res_sigma2c.setConstant(); sig_mean3c.setConstant(); res_sigma3c.setConstant(); frac2c.setConstant(); frac3c.setConstant(); peak_frac1c.setConstant(); //=======================bidimensional fit======================================================== RooProdPdf totalPdf2("totalPdf2","totalPdf",RooArgSet(model3,model3c)) ; TFile *data_file = TFile::Open("Dalitz_Rho_hepg_trig_check_mc_Lcrho.root"); TTree *data_tree = dynamic_cast(data_file->Get("lb;1")); RooDataSet data2("data2","Dataset",data_tree,RooArgSet(mrho46_gRandom,mrho56_gRandom));//,m3pi // data->reduce(Cut("m3pi>2")) ; cout << "=== Data loaded" << endl; RooFitResult *fit= totalPdf2.fitTo(data2,Extended(kFALSE),SumW2Error(kFALSE),Hesse(kTRUE),Save(kTRUE),Minos(kTRUE)); fit->Print("v"); TH2D* pdf_tot = totalPdf2.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 = data2.createHistogram("dh2",mrho46_gRandom,Binning(100),YVar(mrho56_gRandom,Binning(100))); // for data ,Cut("m3pi>2") TCanvas *canvas_plot2Db= new TCanvas("cplot2Db","Plot",1); canvas_plot2Db->Divide(1,2); canvas_plot2Db->cd(1); pdf_tot->Draw("SURFZ"); canvas_plot2Db->cd(2); dh2->Draw("LEGO"); // dh2->Draw("same"); // canvas_plot2D->cd(3); // dh2->Draw("LEGO"); // pdf_tot->Draw("same"); // TCanvas *canvas_fitXb= new TCanvas("cfitxb","Bidimensional fit X Projection",1); RooPlot* framex = mrho46_gRandom.frame() ; data2.plotOn(framex) ; //,Cut("m3pi>2") totalPdf2.plotOn(framex) ; framex->Draw(); TCanvas *canvas_fitYb= new TCanvas("cfityb","Bidimensional fit Y Projection",1); RooPlot* framey = mrho56_gRandom.frame() ; data2.plotOn(framey) ; //,Cut("m3pi>2") totalPdf2.plotOn(framey) ; framey->Draw(); RooWorkspace *ws = new RooWorkspace("Rhopi_ws","Rhopi template WS"); ws->import(totalPdf2,RenameAllVariablesExcept("RhopiVar","mrho46_gRandom,mrho56_gRandom,dm,k"),RenameAllNodes("totalPdf_RhopiPdf")); ws->writeToFile("Rhopi_WS_template2D.root"); }