void fit_toydata() { gROOT->Reset(); gStyle->SetOptFit(kTRUE); Color_t white=10; gStyle->SetCanvasColor(white);gSystem->Load("libRooFit"); using namespace RooFit; using namespace RooStats; std::ostringstream name; std::ostringstream mode; // #ifdef __CINT__ //gROOT->ProcessLineSync(".x nLGn.cxx+"); //#endif RooRealVar *deltam = new RooRealVar("deltam","", 0.139, 0.158); RooRealVar *md0 = new RooRealVar("md0","", 1.8, 1.9); // DeltaM fit RooRealVar *mean_dm = new RooRealVar("mean_dm", "MEAN of 1st gaussian",0.145866,0.14,0.158); RooRealVar *sigma_dm = new RooRealVar("sigma_dm", "Sigma of 1st gaussian",0.001829,0,0.0089); RooGaussian *gauss_dm= new RooGaussian("gauss_dm", "1st gaussian PDF", *deltam, *mean_dm, *sigma_dm); RooRealVar *delm_dm=new RooRealVar("delm_dm", "Difference of Mean",-0.0003872,-0.01,0.01); RooFormulaVar *mean_dm2= new RooFormulaVar("mean_dm2","Mean of 2nd gaussian", "@0+@1",RooArgList(*mean_dm,*delm_dm)); RooRealVar *mean_dm3 = new RooRealVar("mean_dm3", "MEAN of 1st gaussian",0.1454025,0.14,0.158); RooRealVar *sigma_dm2 = new RooRealVar("sigma_dm2", "Sigma of 2nd gaussian",0.000921,0,0.0089); RooRealVar *sigma_dm3 = new RooRealVar("sigma_dm3", "Sigma of 2nd gaussian",0.0002849,0,0.0089); RooGaussian *gauss_dm2= new RooGaussian("gauss_dm2", "1st gaussian PDF", *deltam, *mean_dm3, *sigma_dm2); RooGaussian *gauss_dm3= new RooGaussian("gauss_dm3", "3rd gaussian PDF", *deltam, *mean_dm3, *sigma_dm3); RooRealVar *mean_dm4 = new RooRealVar("mean_dm4", "MEAN of 1st gaussian",0.1454468,0.14,0.158); RooRealVar *sigma_dm4 = new RooRealVar("sigma_dm4", "Sigma of 2nd gaussian",0.0004967,0,0.0089); RooGaussian *gauss_dm4 = new RooGaussian("gauss_dm4", "3rd gaussian PDF", *deltam, *mean_dm4, *sigma_dm4); RooRealVar *a2a_dm1=new RooRealVar("a2a_dm1", " Area2/Area1",0.1990,0.01,1); RooRealVar *a2a_dm2=new RooRealVar("a2a_dm2", " Area2/Area1",0.4540,0.01,1); RooRealVar *sls1_dm=new RooRealVar("sls1_dm","Ratio of Sigma",0.1,15); RooRealVar *srs1_dm=new RooRealVar("srs1_dm","Ratio of Sigma",0.1,10); RooRealVar *a2a_dm=new RooRealVar("a2a_dm", " Area2/Area1",0.556,0.01,1); RooAddPdf *dG_dm3= new RooAddPdf("dG_dm3", " 1st GAuss + 2nd Gauss", RooArgList(*gauss_dm,*gauss_dm2), RooArgList(*a2a_dm1)); RooAddPdf *dG_dm2= new RooAddPdf("dG_dm2", " 1st GAuss + 2nd Gauss", RooArgList(*dG_dm3,*gauss_dm3), RooArgList(*a2a_dm)); RooAddPdf *dG_dm4= new RooAddPdf("dG_dm4", " 1st GAuss + 2nd Gauss", RooArgList(*dG_dm2,*gauss_dm4), RooArgList(*a2a_dm2)); RooRealVar *mean_dmp = new RooRealVar("mean_dmp", "MEAN of 1st gaussian",0.145438); RooRealVar *sigma_dmp = new RooRealVar("sigma_dmp", "Sigma of 1st gaussian",0.000461); RooGaussian *gauss_dmp= new RooGaussian("gauss_dmp", "1st gaussian PDF", *deltam, *mean_dmp, *sigma_dmp); //-------------------------combinatorial PDF in deltaM----------------- RooRealVar th("th","a",0.139); RooRealVar a("a","a", -278.6,-5000,5000); RooRealVar b("b","",-1095,-100000,100000); RooRealVar c("c","",0); RooRealVar d("d","",0); RooGenericPdf *t_bkg = new RooGenericPdf("t_bkg", "(deltam>=th)*pow((deltam-th),3)*exp(a*(deltam-th)+ b*pow((deltam-th),2))", RooArgSet(*deltam,th,a,b)); //-------------------random slow pion----------------------- RooRealVar th1("th1","th1",0.139); RooRealVar a1("a1","a1",-422.64,-2000,2000); RooRealVar b1("b1","b1", 9894,-10000,10000); RooRealVar c1("c1","",0); RooRealVar d1("d1","",0); RooGenericPdf *t_bkg1 = new RooGenericPdf("t_bkg1", "(deltam>=th1)*pow((deltam-th1),2)*exp(a1*(deltam-th1)+ b1*pow((deltam-th1),2))", RooArgSet(*deltam,th1,a1,b1)); RooRealVar * slope1= new RooRealVar("slope1", "Slope of Polynomial",0.33650,-1,1); RooRealVar * slope2= new RooRealVar("slope2", " Slope of Cheb2",-0.21539,-1,1); RooChebychev *chebpol = new RooChebychev("chebpol","Chebshev Polynomial ", *deltam, RooArgList(*slope1,*slope2)); RooRealVar *a2a_dmc=new RooRealVar("a2a_dmc", " Area2/Area1",0.1,1); RooAddPdf *dG_dmc= new RooAddPdf("dG_dmc", " 1st Bkg + 2nd Bkg", RooArgList(*t_bkg,*chebpol), RooArgList(*a2a_dmc)); // M_D0 fit RooRealVar *mean_md = new RooRealVar("mean_md", "MEAN of 1st gaussian",1.864,1.82,1.9); RooRealVar *sigma_md = new RooRealVar("sigma_md", "Sigma of 1st gaussian",0.00652,0,0.0089); RooGaussian *gauss_md= new RooGaussian("gauss_md", "1st gaussian PDF", *md0, *mean_md, *sigma_md); RooRealVar *mean_md1 = new RooRealVar("mean_md1", "MEAN of 1st gaussian",1.864,1.82,1.9); RooRealVar *sigma_md1 = new RooRealVar("sigma_md1", "Sigma of 2nd gaussian",0.004118,0,0.0089); RooGaussian *gauss_md1= new RooGaussian("gauss_md1", "2nd gaussian PDF", *md0, *mean_md1, *sigma_md1); RooRealVar *sls1_md=new RooRealVar("sls1_md","Ratio of Sigma",0.1,20); RooRealVar *srs1_md=new RooRealVar("srs1_md","Ratio of Sigma",0.1,20); RooRealVar *a2a_md=new RooRealVar("a2a_md", " Area2/Area1",0.1902,0.01,1); RooRealVar *a2a_md1=new RooRealVar("a2a_md1", " Area2/Area1",0.438,0.01,1); RooRealVar *delm_md=new RooRealVar("delm_md", "Difference of Mean",0.00037,-0.01,0.01); RooFormulaVar *mean_md2= new RooFormulaVar("mean_md2","Mean of 2nd gaussian", "@0+@1",RooArgList(*mean_md,*delm_md)); RooFormulaVar *sigmaL_md= new RooFormulaVar("sigmaL_md", "Sigma of 2nd gaussian", " @0*@1", RooArgList(*sigma_md,*sls1_md)); RooFormulaVar *sigmaR_md= new RooFormulaVar("sigmaR_md", "Sigma of 2nd gaussian", " @0*@1", RooArgList(*sigma_md,*srs1_md)); RooBifurGauss *gauss_md3= new RooBifurGauss("gauss_md3", "3nd gaussian PDF", *md0, *mean_md2, *sigmaL_md,*sigmaR_md); RooAddPdf *dG_md1= new RooAddPdf("dG_md1", " 1st GAuss + 3rd Gauss", RooArgList(*gauss_md3,*gauss_md), RooArgList(*a2a_md)); RooAddPdf *dG_md2= new RooAddPdf("dG_md2", " 1st GAuss + 3rd Gauss", RooArgList(*dG_md1,*gauss_md1), RooArgList(*a2a_md1)); //-------------------------random slow pion in M_D----------------------- RooRealVar *mean_mdp = new RooRealVar("mean_mdp", "MEAN of 1st gaussian",1.865032); RooRealVar *sigma_mdp = new RooRealVar("sigma_mdp", "Sigma of 1st gaussian",0.004928); RooGaussian *gauss_mdp= new RooGaussian("gauss_mdp", "1st gaussian PDF", *md0, *mean_mdp, *sigma_mdp); //-------------------------combinatorial PDF in M_D----------------- RooRealVar * slope4= new RooRealVar("slope4", "Slope of Polynomial",-0.26226,-1,1); RooRealVar * slope5= new RooRealVar("slope5", " Slope of Cheb2",0.0457,-1,1); RooChebychev *chebpol1 = new RooChebychev("chebpol1","Chebshev Polynomial ", *md0, RooArgList(*slope4,*slope5)); //----------------------------3 body decay in M_D------------------------- RooRealVar * slope6= new RooRealVar("slope6", "Slope of Polynomial",-0.069,-1,1); RooRealVar * slope7= new RooRealVar("slope7", " Slope of Cheb2",-0.6321,-1,1); RooChebychev *chebpol2 = new RooChebychev("chebpol2","Chebshev Polynomial ", *md0, RooArgList(*slope6,*slope7)); RooProdPdf *pdf_signal_signal = new RooProdPdf("pdf_signal_signal","pdf_signal",RooArgList(*dG_dm4,*dG_md2));//2DsignalMC RooProdPdf *pdf_signal = new RooProdPdf("pdf_signal","pdf_signal",RooArgList(*dG_dmc,*chebpol1));//2Dcomb //RooProdPdf *pdf_signal_peaking = new RooProdPdf("pdf_signal_peaking","pdf_signal_peaking",RooArgList(*gauss_dmp,*gauss_mdp));//2Dpeaking RooProdPdf *pdf_slow_pion1 = new RooProdPdf("pdf_slow_pion1","pdf_slow_pion1",RooArgList(*t_bkg1,*gauss_mdp));//random slow pion RooProdPdf *pdf_slow_pion2 = new RooProdPdf("pdf_slow_pion2","pdf_slow_pion2",RooArgList(*gauss_dmp,*chebpol2));//3 body decay RooRealVar *BKG = new RooRealVar("BKG", "",500,1000000); RooRealVar *bkg_peaking = new RooRealVar("bkg_peaking", "",1465); RooRealVar *sig = new RooRealVar("sig", "",5000,3000000); RooRealVar *spion1 = new RooRealVar("spion1", "",5555,1,500000);//5555); RooRealVar *spion2 = new RooRealVar("spion2", "",3866,1,100000); RooAddPdf *total = new RooAddPdf("total","total",RooArgList(*pdf_signal_signal,*pdf_signal,*pdf_slow_pion1,*pdf_slow_pion2),RooArgList(*sig,*BKG,*spion1,*spion2)); float tempe,temps,tempds,tempde,temp2,temp3,temp4,tempdeb,tempdsb,tempeb,tempsb,tempdl,tempdh ; int Enty = 0; int Entx = 0; TCanvas c3("c3", "First canvas", 800, 800); c3.Divide(2,2); TCanvas c2("c2", "Second canvas", 800, 800); c2.Divide(2,2); TPostScript ps ("toy_pmu-.ps",112); // TPostScript ps1 ("md0_toy_pmu-.ps",112); FILE *f1, *f2; f1 = fopen("bkg-toy.txt","w"); f2 = fopen("signal-toy.txt","w"); for ( int i=1; i<2; i++){ string infiles="./incbkg"; std::ostringstream temp1; temp1 << infiles << "/bkginc-" << i; infiles = temp1.str(); infiles.append(".txt"); std::cout << infiles << endl; //RooDataSet *data = RooDataSet::read(infiles.c_str(), RooArgList(*deltam,*md0)); RooDataHist *data = RooDataHist::read(infiles.c_str(), RooArgList(*deltam,*md0)); string infile1="./incdat"; std::ostringstream temp2; temp2 << infile1 << "/siginc-" << i; infile1 = temp2.str(); infile1.append(".txt"); std::cout << infile1 << endl; //RooDataSet *data1 = RooDataSet::read(infile1.c_str(), RooArgList(*deltam,*md0)); RooDataHist *data1 = RooDataHist::read(infile1.c_str(), RooArgList(*deltam,*md0)); string infile2="./spion1"; std::ostringstream temp3; temp3 << infile2 << "/spion1-" << i; infile2 = temp3.str(); infile2.append(".txt"); std::cout << infile2 << endl; //RooDataSet *data2 = RooDataSet::read(infile2.c_str(), RooArgList(*deltam,*md0)); RooDataHist *data2 = RooDataHist::read(infile2.c_str(), RooArgList(*deltam,*md0)); string infile3="./spion2"; std::ostringstream temp4; temp4 << infile3 << "/spion2-" << i; infile3 = temp4.str(); infile3.append(".txt"); std::cout << infile3 << endl; //RooDataSet *data3 = RooDataSet::read(infile3.c_str(), RooArgList(*deltam,*md0)); RooDataHist *data3 = RooDataHist::read(infile3.c_str(), RooArgList(*deltam,*md0)); data->append(*data1); data->append(*data2); data->append(*data3); RooFitResult *fitRes = total->fitTo(*data,Extended(true),Minos(true),Save()); RooPlot* xframe = deltam->frame(Bins(40)); data->plotOn(xframe); // total->paramOn(xframe); total->plotOn(xframe); //dG->plotOn(deltaEplot,Components(RooArgSet(*dG_2)),LineColor(kGreen),LineStyle(kDashed)); //dG->plotOn(deltaEplot,Components(RooArgSet(*lg)),LineColor(kRed),LineStyle(kDashed)); total->plotOn(xframe,Components(RooArgSet(*total)),LineColor(kGreen),LineStyle(kDashed)); name.str(""); name << "#DeltaM [" << i <<"]"; xframe->SetTitle(name.str().c_str()); //c3->SaveAs(finaF.c_str()); if (Entx==0) {ps.NewPage();} // if (Entx==3 ) {c2.cd(4); Entx=0; //xframe->Draw(); // c2.Update();} //else { c2.cd(Entx+1); std::cout<<"Entx"<Draw(); Entx++; //} tempds = sig->getVal(); tempde = sig->getError(); tempdl = sig->getAsymErrorLo(); tempdh = sig->getAsymErrorHi(); fprintf(f2,"%f \t %f \t %f \t %f\n",tempds,tempde,tempdl,tempdh); RooPlot* yframe = md0->frame(Bins(40)); data->plotOn(yframe); total->paramOn(yframe); total->plotOn(yframe); //dG->plotOn(deltaEplot,Components(RooArgSet(*dG_2)),LineColor(kGreen),LineStyle(kDashed)); //dG->plotOn(deltaEplot,Components(RooArgSet(*lg)),LineColor(kRed),LineStyle(kDashed)); total->plotOn(yframe,Components(RooArgSet(*total)),LineColor(kGreen),LineStyle(kDashed)); name.str(""); name << "M_{D^{0}} [" << i <<"]"; yframe->SetTitle(name.str().c_str()); //c3->SaveAs(finaF.c_str()); //if (Enty==0) {ps1.NewPage();} if (Entx==3 ) {c2.cd(4); Entx=0; yframe->Draw(); c2.Update();} else {c2.cd(Entx+1); Entx++; yframe->Draw();} tempsb = BKG->getVal(); tempeb = BKG->getError(); fprintf(f1,"%f \t %f \n",tempsb,tempeb); //fprintf(f2,"%f \t %f \n",temp3,temp4); } fclose (f1); fclose (f2); ps.Close(); // ps1.Close(); return 0; }