using namespace RooFit ; using namespace std ; int method1() { RooMsgService::instance().getStream(0).removeTopic(Eval); RooMsgService::instance().getStream(1).removeTopic(Eval); //---------------declare variables for Mds pdf signal and background shape parameters-----// gROOT->LoadMacro("/Users/amansangal/Downloads/recon2_newbcs2/style_file/style.C"); //-----------------3 fit variables------------// RooRealVar t("t","t_{rec} (ns)",-0.002,0.004); //decay time t RooRealVar dt("dt","#sigma_{t} (ns)",0,0.0009); // error in decay time t RooDataSet* data1 = (RooDataSet*)RooDataSet::read("signal_2dSR_1.txt",RooArgSet(t,dt),"Q"); RooDataSet* data2 = (RooDataSet*)RooDataSet::read("signal_2dSR_2.txt",RooArgSet(t,dt),"Q"); data1->append(*data2); data1->Print("V"); //--------frame for decay time t-------------// RooPlot* framet = t.frame(); data1->plotOn(framet,Binning(100)); //-------------frame for error in decay time dt------------------------------------// RooPlot* framedt = dt.frame(); data1->plotOn(framedt,Binning(80)); //----------------------BLOCK 4------------------------// //--------here start constructing the models to be used for fitting -----------// //--------Part I Signal-----------// //----------Resolution function for signal--------------// //--------------gauss1---------------// RooRealVar mean("mean","mean",0.0001,-2,0.007); // mean/bias RooRealVar sc1("sc1","sc1",1.0,0,2); // scaling factor // RooFormulaVar sigma_t("sigma_t","sc1*dt",RooArgList(sc1,dt)); // RooGaussModel gauss1("gauss1","gauss",t,mean,sigma_t); // include decay time error RooGaussModel gauss1("gauss1","gauss",t,mean,sc1,dt); //----------------gauss2-----------// RooRealVar sc2("sc2","sc2",1.9933,0,10); // scaling factor // RooFormulaVar sigma_t1("sigma_t1","sc2*dt",RooArgList(sc2,dt)); // RooGaussModel gauss2("gauss2","gauss",t,mean,sigma_t1); // include decay time error RooGaussModel gauss2("gauss2","gauss",t,mean,sc2,dt); //------------gauss1+gauss2---------// RooRealVar fg1r("fg1r","fg1r",0.965191,0,1); RooAddModel g2("g2","g2",RooArgList(gauss1,gauss2),fg1r); //-------lifetime variable---------// RooRealVar tau("tau","tau",0.000104,0,0.001); //-------------final convolved pdf for signal--------// RooDecay decay("decay","lifetime",t,tau,gauss1,RooDecay::SingleSided); // only 1G as resolution function //--------------pdf for dt signal--------// //--------------Johnson Su function -----------------// RooRealVar xi("xi","xi",0.00002,0,0.00006); RooRealVar lambda("lambda","lambda",0.00003,0,0.0003); RooRealVar delta("delta","delta",1,0,5); RooRealVar gamma("gamma","gamma",-2,-20,2); RooJohnson jsu("jsu","jsu",dt,xi,lambda,gamma,delta); //---------gaussian---------// RooRealVar mean_dt("mean_dt","mean",0.0001,0,0.0006); RooRealVar sigma_dt("sigma_dt","sigma",0.00004,0,0.0008); RooGaussian gauss1_dt("gauss1_dt","gauss1",dt,mean_dt,sigma_dt); RooRealVar fr_g_dt("fr_g_dt","fr_g_dt",0.5,0,1); RooAddPdf gjsu("gjsu","gjsu",RooArgList(gauss1_dt,jsu),RooArgList(fr_g_dt)); //------final signal pdf --------------// RooProdPdf sig_3d("sig_3d","sig_3d",RooArgSet(gjsu),Conditional(decay,t)); sig_3d.fitTo(*data1,Minos(0),Save(),Timer(kTRUE),Hesse(0),NumCPU(8)); RooFitResult* fitres = sig_3d.fitTo(*data1,Minos(1),Save(),Timer(kTRUE),Hesse(0),NumCPU(8)); sig_3d.plotOn(framet,LineColor(kBlue),ProjWData(dt,*data1)); // finalp.plotOn(framet,LineColor(kBlue)); sig_3d.paramOn(framet); // finalp.plotOn(framet,Components(RooArgList(bkg_3d)),LineColor(kRed),LineStyle(7), LineWidth(2)); sig_3d.plotOn(framedt,LineColor(kBlue)); // finalp.plotOn(framedt,Components(RooArgList(bkg_3d)),LineColor(kRed),LineStyle(7), LineWidth(2)); //--------------getting # of parameters-------------// RooArgSet observables(t,dt); RooArgSet *flparams = sig_3d.getParameters(observables); int nparam = (flparams->selectByAttrib("Constant",kFALSE))->getSize(); cout<<" ndf are "<chiSquare(nparam);// cout<<"chi square = "<SetFillColor(10); box->SetBorderSize(1); box->SetTextAlign(12); box->SetTextSize(0.06F); 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 ); framet->addObject(box) ; /////////////pull distribution//////////////// //****************************** RooPlot* z1frame = t.frame(Title("Pull Distribution")); RooHist* hpull1 = framet->pullHist(); hpull1->SetFillColor(4); hpull1->SetLineColor(0); z1frame->addPlotable( hpull1, "B" ); TCanvas* c1 = new TCanvas("c1","c1",550,400) ; c1->Divide(1,2) ; // column, row double xmin = -0.002; double xmax = 0.004; 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); c1->cd(2)->SetPad(0.005,0.005,0.995,0.2525); line->SetLineColor(kRed); line->SetLineWidth(2); z1frame->GetYaxis()->SetTitle("Pull"); z1frame->GetYaxis()->SetTitleSize(0.18); z1frame->GetYaxis()->SetTitleOffset(0.37); z1frame->GetYaxis()->SetLabelSize(0.08); z1frame->GetXaxis()->CenterTitle(); z1frame->Draw(); line1->SetLineColor(kRed); line1->SetLineWidth(3); line1->SetLineStyle(2); line2->SetLineColor(kRed); line2->SetLineWidth(3); line2->SetLineStyle(2); line->Draw("SAME"); //z1frame->GetYaxis()->SetRangeUser(-3.5, 3.5); line1->Draw("SAME"); line2->Draw("SAME"); c1->cd(1)->SetPad(0.005,0.2525,0.995,0.995); //gPad->SetLeftMargin(0.15); framet->GetXaxis()->CenterTitle(); // sig_3d.plotOn(framet,Components(RooArgList(bkg_3d)),LineColor(kRed),LineStyle(7), LineWidth(2)); //sig_3d.plotOn(framet,Components(RooArgList(gaussmb)),LineColor(kRed),LineStyle(7), LineWidth(2)); framet->Draw(); c1->Draw(); c1->SaveAs("tproj.C"); c1->SaveAs("tproj.eps"); TCanvas *cdt = new TCanvas("cdt","cdt",550,400); cdt->Divide(1,2) ; //--------------dt pull plot---------// double xmindt = 0; double xmaxdt = 0.0009; TLine *linedt = new TLine(xmindt,0.0,xmaxdt,0.0); TLine *line1dt = new TLine(xmindt,-3.0,xmaxdt,-3.0); TLine *line2dt = new TLine(xmindt,3.0,xmaxdt,3.0); RooPlot* z2frame = dt.frame(Title("Pull Distribution")); RooHist* hpull2 = framedt->pullHist(); hpull2->SetFillColor(4); hpull2->SetLineColor(0); z2frame->addPlotable( hpull2, "B" ); cdt->cd(2)->SetPad(0.005,0.005,0.995,0.2525); line->SetLineColor(kRed); line->SetLineWidth(2); z2frame->GetYaxis()->SetTitle("Pull"); z2frame->GetYaxis()->SetTitleSize(0.18); z2frame->GetYaxis()->SetTitleOffset(0.37); z2frame->GetYaxis()->SetLabelSize(0.08); z2frame->GetXaxis()->CenterTitle(); z2frame->Draw(); line1dt->SetLineColor(kRed); line1dt->SetLineWidth(3); line1dt->SetLineStyle(2); line2dt->SetLineColor(kRed); line2dt->SetLineWidth(3); line2dt->SetLineStyle(2); linedt->Draw("SAME"); //z2frame->GetYaxis()->SetRangeUser(-3.5, 3.5); line1dt->Draw("SAME"); line2dt->Draw("SAME"); cdt->cd(1)->SetPad(0.005,0.2525,0.995,0.995); framedt->Draw(); cdt->Draw(); cdt->SaveAs("dtproj.C"); cdt->SaveAs("dtproj.eps"); fitres->Print("v"); flparams->printLatex() ; flparams->printLatex(OutputFile("fitresult.tex")) ; ofstream param_file; param_file.open ("example.txt"); // delta,fr_g_t,gamma,lambda,mean,mean_dt,sc1,sigma_dt,tau,xi param_file<<"Parameter Name "<<" Fit Result"<Write(); cdt->Write(); fcan->Close(); return 1; }