{ using namespace RooFit; using namespace RooStats; gStyle->SetOptFit(kTRUE); Color_t white=10; gStyle->SetCanvasColor(white); gSystem->Load("libRooFit"); RooRealVar *delte=new RooRealVar("delte","delte",-0.09,0.09); RooRealVar *ymass=new RooRealVar("ymass","ymass",3.835,3.91); RooRealVar *mbc=new RooRealVar("mbc","mbc",5.27,5.29); string finaO="del.txt"; string finaOm="mass.txt"; fstream likh(finaO.c_str(), ios::out| ios::app); fstream likhm(finaOm.c_str(), ios::out| ios::app); //Loop for(int i=1;i<2;i++) { cout<<":::::::::::::::::::::::::::::::::::START::::::::::::::::::::::::::::::::::::::::"<setVal(x_delte); ymass->setVal(x_ymass); mbc->setVal(x_mbc); data->add(RooArgList(*delte,*ymass,*mbc)); } //___________________________________________________________ //DeltaE PDF //_________________________________________________________ RooRealVar * slope1= new RooRealVar("slope1", "Slope of Polynomial", -0.1628,-11,11); RooRealVar * slope2= new RooRealVar("slope2", "Slope of Polynomial", 0.3807, -11,11); RooChebychev *chebpol = new RooChebychev("chebpol","Chebshev Polynomial ", *delte, RooArgList(*slope1,*slope2)); RooRealVar *mean = new RooRealVar("mean", "MEAN of 1st gaussian",-0.0005104, -0.01,0.01); RooRealVar *sigma1=new RooRealVar("sigma1", "Sigma of 1st gaussian", 0.00736,0.00001,0.02); RooRealVar *delm=new RooRealVar("delm", "Difference of Mean",0); RooRealVar *s2s1=new RooRealVar("s2s1","Ratio of Sigma",0.629); RooFormulaVar *mean2=new RooFormulaVar("mean3","Mean of second gaussian","@0+@1",RooArgList(*mean,*delm)); RooFormulaVar *sigma2=new RooFormulaVar("sigma3","sigma","@0*@1",RooArgList(*sigma1,*s2s1)); RooRealVar *a3a=new RooRealVar("a3a", " Area2/Area1",0.1363); RooRealVar *a2a=new RooRealVar("a2a", " Area2/Area1",0.554); RooRealVar *delm1=new RooRealVar("delm1", "Difference of Mean",0.00341); RooRealVar *sls1=new RooRealVar("sls1","Ratio of Sigma", 3.105); RooRealVar *srs1=new RooRealVar("srs1","Ratio of Sigma", 2.199); RooFormulaVar *mean3= new RooFormulaVar("mean2","Mean of 2nd gaussian", "@0+@1",RooArgList(*mean,*delm1)); RooFormulaVar *sigmaL= new RooFormulaVar("sigmaL", "Sigma of 2nd gaussian", " @0*@1", RooArgList(*sigma1,*sls1)); RooFormulaVar *sigmaR= new RooFormulaVar("sigmaR", "Sigma of 2nd gaussian", " @0*@1", RooArgList(*sigma1,*srs1)); RooGaussian *gauss1_1= new RooGaussian("gauss1_1", "1st gaussian PDF", *delte, *mean, *sigma1); RooGaussian *gauss2_1= new RooGaussian("gauss2_1", "3rd gaussian PDF", *delte, *mean2, *sigma2); RooBifurGauss *gauss3_1= new RooBifurGauss("gauss3_1", "2nd gaussian PDF", *delte, *mean3, *sigmaL,*sigmaR); RooAddPdf *dG= new RooAddPdf("dG", " 1st GAuss + 2nd Gauss", RooArgList(*gauss3_1,*gauss2_1,*gauss1_1), RooArgList(*a3a,*a2a)); RooRealVar *sig = new RooRealVar("sig", "",190,10,2000); RooRealVar *bkg = new RooRealVar("bkg", "",100,0,5000); RooAddPdf *depdf = new RooAddPdf ("depdf", "Two Gaussian + ",RooArgList(*dG,*chebpol), RooArgList(*sig,*bkg)); //_________________________________________________ //Mass of B(mbc) //_________________________________________________ RooRealVar *meanm=new RooRealVar("meanm","Mean of ist gaussian",5.278999,5.275,5.279); RooRealVar *sigmam=new RooRealVar("sigmam","sigma",0.0025982,0.001,0.004); RooGaussian *mgau=new RooGaussian("mgau"," ",*mbc,*meanm,*sigmam); RooRealVar *Margpar=new RooRealVar("Margpar","argus shape",-224,-300,-1); RooArgusBG *argus=new RooArgusBG("argus","ArgusPDF",*mbc, RooConst(5.289),*Margpar); RooProdPdf *restbkg=new RooProdPdf("restbkg","",RooArgSet(*chebpol,*argus)); RooProdPdf *msig= new RooProdPdf ("msig","",RooArgSet(*dG,*mgau)); RooAddPdf *tot=new RooAddPdf("tot","",RooArgSet(*msig,*restbkg),RooArgList(*sig,*bkg)); RooRealVar *sope1=new RooRealVar("sope1","slope",-10,10); RooRealVar *sope2=new RooRealVar("sope2","slope",-10,10); RooChebychev *mpol=new RooChebychev("mpol","chebychev",*mbc,RooArgList(*sope1,*sope2)); RooRealVar *mbsig=new RooRealVar("mbsig","",785100,1000,400000); RooRealVar *mBKG=new RooRealVar("mBKG","",7851000,-100,40000000); RooAddPdf *mpdf=new RooAddPdf("mpdf","Gaussian+",RooArgList(*mgau,*mpol),RooArgList(*mbsig,*mBKG)); //_________________________________________________ //Mass //_________________________________________________ RooRealVar * xope1= new RooRealVar("xope1", "Slope of Polynomial", -0.4704,-5,5); RooRealVar * xope2= new RooRealVar("xope2", "Slope of Polynomial", -0.3888,-11,11); RooChebychev *xpol = new RooChebychev("xpol","Chebshev Polynomial ", *ymass, RooArgList(*xope1,*xope2)); RooRealVar *mean_X=new RooRealVar("mean_X","Mean of first gauusian",3.8713,3.835,3.9); RooRealVar *sigma1_X=new RooRealVar("sigma1_X","Sigma of first gaussain",0.003679,0,0.3); RooRealVar *a2a_X=new RooRealVar("a2a_X","Area2/Area1", 0.598); RooRealVar *delm_X=new RooRealVar("delm_X","difference of means",0); RooRealVar *s2s1_X=new RooRealVar("s2s1_X","Ratio of sigma",0.4733); RooFormulaVar *mean2_X=new RooFormulaVar("mean2_X","Mean of second gaussian","@0+@1", RooArgList(*mean_X,*delm_X)); RooFormulaVar *sigma2_X=new RooFormulaVar("sigma2_X","","@0*@1", RooArgList(*sigma1_X,*s2s1_X)); RooRealVar *a3a_X=new RooRealVar("a3a","Area3/Area1",0.0949); RooRealVar *delm1_X=new RooRealVar("delm1","difference of means",-0.00084); RooRealVar *sls1_X=new RooRealVar("sls1","Ratio of sigma",3.059); RooRealVar *srs1_X=new RooRealVar("srs1","Ratio of sigma",3.028); RooFormulaVar *mean3_X=new RooFormulaVar("mean3","Mean of second gaussian","@0+@1", RooArgList(*mean_X,*delm1_X)); RooFormulaVar *sigmaL_X=new RooFormulaVar("sigmaL_X","","@0*@1", RooArgList(*sigma1_X,*sls1_X)); RooFormulaVar *sigmaR_X=new RooFormulaVar("sigmaR_X","","@0*@1", RooArgList(*sigma1_X,*srs1_X)); RooGaussian *gauss1_X= new RooGaussian("gauss1_X", "1st gaussian PDF", *ymass, *mean_X, *sigma1_X); RooGaussian *gauss2_X= new RooGaussian("gauss2_X", "2nd gaussian PDF", *ymass, *mean2_X, *sigma2_X); RooBifurGauss *gauss3_X= new RooBifurGauss("gauss3_X", "3rd gaussian PDF", *ymass, *mean3_X, *sigmaL_X,*sigmaR_X); RooAddPdf *dG_X= new RooAddPdf("dG_X", " 1st GAuss + 2nd Gauss", RooArgList(*gauss3_X,*gauss2_X, *gauss1_X), RooArgList(*a3a_X,*a2a_X)); RooRealVar *xsig = new RooRealVar("xsig", "",187,25,500); RooRealVar *xBKG = new RooRealVar("xBKG", "",100,-10,2000); RooAddPdf *ypdf = new RooAddPdf ("ypdf", "Two Gaussian + ",RooArgList(*dG_X,*xpol), RooArgList(*xsig,*xBKG)); RooFitResult *fitRes = depdf->fitTo(*data,Extended(true),Minos(true),Save()); string finaF; std::ostringstream teF; teF << "./fig/Sig-_" << i; finaF = teF.str(); finaF.append(".gif"); TCanvas* c1 = new TCanvas("c1", "", 900, 700); gStyle->SetOptStat(0); gPad -> SetFillColor(0); RooPlot* deltaEplot = delte->frame(Bins(50)); data->plotOn(deltaEplot); depdf->paramOn(deltaEplot); depdf->plotOn(deltaEplot); depdf->plotOn(deltaEplot,Components(RooArgSet(*chebpol)),LineColor(kGreen),LineStyle(kDashed)); depdf->plotOn(deltaEplot,Components(RooArgSet(*dG)),LineColor(kRed),LineStyle(kDashed)); deltaEplot->Draw(); c1->SaveAs(finaF.c_str()); likh << i << "\t" << fitRes->status() << "\t" << sig->getVal() << "\t" <getError() <<"\t" << sig->getAsymErrorHi() << "\t" << sig->getAsymErrorLo() << "\t" << bkg->getVal() << "\t" << bkg->getError() << "\t" << bkg->getAsymErrorHi() << "\t" << bkg->getAsymErrorLo() << endl; cout<<":::::::::::::::::::::::::::::::::::END DeltaE::::::::::::::::::::::::::::::::::::::::"<setConstant(true); sigma1->setConstant(true); slope1->setConstant(true); slope2->setConstant(true); bkg->setConstant(true); RooStats::SPlot* sData = new RooStats::SPlot("sWeightedData","sWeightedData", *data, tot, RooArgList(*sig,*bkg) ); std::cout << " Check Weights \n"; std::cout << "\n Yield of signale is " << sig->getVal() << ". From sWeight it is " << sData->GetYieldFromSWeight("sig") << endl; std::cout << "\n Yield of BKG is " << bkg->getVal() << ". From sWeight it is " << sData->GetYieldFromSWeight("bkg") << endl; RooRealVar *weight = new RooRealVar("weight","my weight", -2000, 20000); RooDataSet *dataS = new RooDataSet("dataS", "data",RooArgSet(*delte,*ymass, *mbc,*weight)); RooDataSet *iData = new RooDataSet("iData", "data",RooArgSet(*ymass,*weight)); RooDataSet *iData0 = new RooDataSet("iData0", "data",RooArgSet(*ymass,*weight)); std::cout << " HEREEE" << dataS->numEntries() ; dataS->Print(); for(int k = 0; k < data->numEntries(); k++){ double val = data->get(k)->getRealValue("delte"); double val2 = data->get(k)->getRealValue("mbc"); double val3 = data->get(k)->getRealValue("ymass"); double we = sData->GetSWeight(k,"sig"); delte ->setVal( val ); mbc ->setVal(val2); ymass->setVal(val3); weight->setVal(we); dataS->add( RooArgList(*delte, *ymass, *mbc, *weight)); // if (we >0) iData->add(RooArgList(*ymass,*weight)); if (we >0) iData0->add(RooArgList(*ymass,*weight)); } TCanvas *c6 =new TCanvas ("c6", "mass sPlot", 900, 600); dataS->Print(); RooFormulaVar wFunc6("wFunc6","event weight","@1",RooArgSet(*ymass,*weight)); RooRealVar* w6 = (RooRealVar*) dataS->addColumn(wFunc6); RooDataSet* wdata6 = new RooDataSet(dataS->GetName(),dataS->GetTitle(),dataS,*dataS->get(),0,w6->GetName()); std::cout << "PRINT " << endl; wdata6->Print(); c6->Draw(); c6->cd(); RooPlot * frame6 = ymass->frame(Bins(100)); wdata6->plotOn(frame6, DataError(RooAbsData::SumW2)); frame6->Draw(); RooDataSet *iwdata = new RooDataSet("iwdata", "", *iData->get(), Import(*iData), WeightVar("weight")); TCanvas *c3 =new TCanvas ("c3", "mass sPlot", 700, 600); RooFitResult *fitresult1 = ypdf->fitTo(*wdata6,Save(),Extended(kTRUE), SumW2Error(kFALSE), Minos(true)); RooPlot* xplotbc = ymass->frame(Bins(50)); wdata6->plotOn(xplotbc, DataError(RooAbsData::SumW2)); ypdf->paramOn(xplotbc); ypdf->plotOn(xplotbc); ypdf->plotOn(xplotbc,Components(RooArgSet(*xpol)),LineColor(kGreen), LineStyle(kDashed)); ypdf->plotOn(xplotbc,Components(RooArgSet(*dG_X)),LineColor(kRed), LineStyle(kDashed)); xplotbc->Draw(); likhm << i << "\t" << fitresult1->status() << "\t" << xsig->getVal() << "\t" <getError() <<"\t" << xsig->getAsymErrorHi() << "\t" << xsig->getAsymErrorLo() << "\t" << xBKG->getVal() << "\t" << xBKG->getError() << "\t" << xBKG->getAsymErrorHi() << "\t" << xBKG->getAsymErrorLo() << endl; string mfinaF; std::ostringstream mteF; mteF << "./mfig/Sig-_" << i; mfinaF = mteF.str(); mfinaF.append(".gif"); c3->SaveAs(mfinaF.c_str()); cout<<":::::::::::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::::::::"<