using namespace RooFit; using namespace RooStats; void jan_2spinAsy_constrain(double nObs1=2, double nObs2=4, double Amax=0.6){ RooWorkspace* ws = new RooWorkspace("jan2Asy"); // variables: predicetd yields RooRealVar* n1 = (RooRealVar*) ws->factory("n1[15,0,50]"); RooRealVar* n2 = (RooRealVar*) ws->factory("n2[15,0,50]"); // model parameterization RooRealVar* A= (RooRealVar*) ws->factory(Form("A[0.12,%.2f,%.2f]",-Amax,Amax)); RooRealVar* N0= (RooRealVar*) ws->factory("N0[15,0,50]"); /* decide what role plays each variable */ ws->defineSet("obs","n1,n2"); ws->defineSet("parOfInterest", "A"); ws->defineSet("nuisPar", "N0"); // intermediate variables depend on the true params of the model ws->factory("expr::mu1('N0*(1+A)',N0,A)"); ws->factory("expr::mu2('N0*(1-A)',N0,A)"); // define pdf's RooAbsPdf* model1 = (RooAbsPdf*) ws->factory("Poisson::pdf1(n1,mu1)"); RooAbsPdf* model2 = (RooAbsPdf*) ws->factory("Poisson::pdf2(n2,mu2)"); ws->factory("PROD::pdfTotal(pdf1,pdf2)"); // First create empty data set, set data value and put into workspace RooDataSet* dataS = new RooDataSet("dataJan", "data 2spin", RooArgSet(*n1,*n2)); // empty // add result of experiment: a pair of 2 measurements for + & - spin state n1->setVal(nObs1); n2->setVal(nObs2); dataS->add(RooArgSet(*n1,*n2)); dataS->Print("v"); ws->import(*dataS); ws->Print(); // return; // Make ModelConfig object and put in workspace ModelConfig* modelConfig = new ModelConfig("recoSSA"); modelConfig->SetWorkspace(*ws); modelConfig->SetPdf("pdfTotal"); modelConfig->SetParametersOfInterest(*ws->set("parOfInterest")); modelConfig->SetNuisanceParameters(*ws->set("nuisPar")); ws->import(*modelConfig); ProfileLikelihoodCalculator plc(*dataS, *modelConfig); double CL=0.6827; plc.SetConfidenceLevel(CL); LikelihoodInterval* plInt = plc.GetInterval(); double Lo = plInt->LowerLimit(*A); double Up = plInt->UpperLimit(*A); TString outTxt=Form("ProfLike CL=%.4f --> A limit[%.3f, %.3f]\n",CL,Lo,Up); cout<Gauss method recA=%.3f 1sigma interval[%.3f, %.3f]\nA\n",nObs1,nObs2,recoA,recoA-asyErr,recoA+asyErr); // return; // PLOT likelihood for 1-sigma c=new TCanvas(); LikelihoodIntervalPlot* intPlot = new LikelihoodIntervalPlot(plInt); intPlot->Draw(); tx=new TText(-0.5,1.6,outTxt); tx->Draw(); }