Zbi_Zgamma.C: Demonstraite Z_Bi = Z_Gamma | Roostats tutorials | rs102_hypotestwithshapes.C: rs102_hypotestwithshapes for RooStats project |
///////////////////////////////////////////////////////////////////////// // // 'Limit Example' RooStats tutorial macro #101 // author: Kyle Cranmer // date June. 2009 // // This tutorial shows an example of creating a simple // model for a number counting experiment with uncertainty // on both the background rate and signal efficeincy. We then // use a Confidence Interval Calculator to set a limit on the signal. // // ///////////////////////////////////////////////////////////////////////// #ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooProfileLL.h" #include "RooAbsPdf.h" #include "RooStats/HypoTestResult.h" #include "RooRealVar.h" #include "RooPlot.h" #include "RooDataSet.h" #include "RooTreeDataStore.h" #include "TTree.h" #include "TCanvas.h" #include "TLine.h" #include "TStopwatch.h" #include "RooStats/ProfileLikelihoodCalculator.h" #include "RooStats/MCMCCalculator.h" #include "RooStats/UniformProposal.h" #include "RooStats/FeldmanCousins.h" #include "RooStats/NumberCountingPdfFactory.h" #include "RooStats/ConfInterval.h" #include "RooStats/PointSetInterval.h" #include "RooStats/LikelihoodInterval.h" #include "RooStats/LikelihoodIntervalPlot.h" #include "RooStats/RooStatsUtils.h" #include "RooStats/ModelConfig.h" #include "RooStats/MCMCInterval.h" #include "RooStats/MCMCIntervalPlot.h" #include "RooStats/ProposalFunction.h" #include "RooStats/ProposalHelper.h" #include "RooFitResult.h" // use this order for safety on library loading using namespace RooFit ; using namespace RooStats ; void rs101_limitexample() { ///////////////////////////////////////// // An example of setting a limit in a number counting experiment with uncertainty on background and signal ///////////////////////////////////////// // to time the macro TStopwatch t; t.Start(); ///////////////////////////////////////// // The Model building stage ///////////////////////////////////////// RooWorkspace* wspace = new RooWorkspace(); wspace->factory("Poisson::countingModel(obs[150,0,300], sum(s[50,0,120]*ratioSigEff[1.,0,2.],b[100,0,300]*ratioBkgEff[1.,0.,2.]))"); // counting model // wspace->factory("Gaussian::sigConstraint(ratioSigEff,1,0.05)"); // 5% signal efficiency uncertainty // wspace->factory("Gaussian::bkgConstraint(ratioBkgEff,1,0.1)"); // 10% background efficiency uncertainty wspace->factory("Gaussian::sigConstraint(1,ratioSigEff,0.05)"); // 5% signal efficiency uncertainty wspace->factory("Gaussian::bkgConstraint(1,ratioBkgEff,0.1)"); // 10% background efficiency uncertainty wspace->factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)"); // product of terms wspace->Print(); RooAbsPdf* modelWithConstraints = wspace->pdf("modelWithConstraints"); // get the model RooRealVar* obs = wspace->var("obs"); // get the observable RooRealVar* s = wspace->var("s"); // get the signal we care about RooRealVar* b = wspace->var("b"); // get the background and set it to a constant. Uncertainty included in ratioBkgEff b->setConstant(); RooRealVar* ratioSigEff = wspace->var("ratioSigEff"); // get uncertaint parameter to constrain RooRealVar* ratioBkgEff = wspace->var("ratioBkgEff"); // get uncertaint parameter to constrain RooArgSet constrainedParams(*ratioSigEff, *ratioBkgEff); // need to constrain these in the fit (should change default behavior) // Create an example dataset with 160 observed events obs->setVal(160.); RooDataSet* data = new RooDataSet("exampleData", "exampleData", RooArgSet(*obs)); data->add(*obs); RooArgSet all(*s, *ratioBkgEff, *ratioSigEff); // not necessary modelWithConstraints->fitTo(*data, RooFit::Constrain(RooArgSet(*ratioSigEff, *ratioBkgEff))); // Now let's make some confidence intervals for s, our parameter of interest RooArgSet paramOfInterest(*s); ModelConfig modelConfig(new RooWorkspace()); modelConfig.SetPdf(*modelWithConstraints); modelConfig.SetParametersOfInterest(paramOfInterest); // First, let's use a Calculator based on the Profile Likelihood Ratio //ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest); ProfileLikelihoodCalculator plc(*data, modelConfig); plc.SetTestSize(.05); ConfInterval* lrint = plc.GetInterval(); // that was easy. // Let's make a plot TCanvas* dataCanvas = new TCanvas("dataCanvas"); dataCanvas->Divide(2,1); dataCanvas->cd(1); LikelihoodIntervalPlot plotInt((LikelihoodInterval*)lrint); plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S"); plotInt.Draw(); // Second, use a Calculator based on the Feldman Cousins technique FeldmanCousins fc(*data, modelConfig); fc.UseAdaptiveSampling(true); fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed fc.SetNBins(100); // number of points to test per parameter fc.SetTestSize(.05); // fc.SaveBeltToFile(true); // optional ConfInterval* fcint = NULL; fcint = fc.GetInterval(); // that was easy. RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true)); // Third, use a Calculator based on Markov Chain monte carlo // Before configuring the calculator, let's make a ProposalFunction // that will achieve a high acceptance rate ProposalHelper ph; ph.SetVariables((RooArgSet&)fit->floatParsFinal()); ph.SetCovMatrix(fit->covarianceMatrix()); ph.SetUpdateProposalParameters(true); ph.SetCacheSize(100); ProposalFunction* pdfProp = ph.GetProposalFunction(); // that was easy MCMCCalculator mc(*data, modelConfig); mc.SetNumIters(20000); // steps to propose in the chain mc.SetTestSize(.05); // 95% CL mc.SetNumBurnInSteps(40); // ignore first N steps in chain as "burn in" mc.SetProposalFunction(*pdfProp); mc.SetLeftSideTailFraction(0.5); // find a "central" interval MCMCInterval* mcInt = (MCMCInterval*)mc.GetInterval(); // that was easy // Get Lower and Upper limits from Profile Calculator cout << "Profile lower limit on s = " << ((LikelihoodInterval*) lrint)->LowerLimit(*s) << endl; cout << "Profile upper limit on s = " << ((LikelihoodInterval*) lrint)->UpperLimit(*s) << endl; // Get Lower and Upper limits from FeldmanCousins with profile construction if (fcint != NULL) { double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s); double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s); cout << "FC lower limit on s = " << fcll << endl; cout << "FC upper limit on s = " << fcul << endl; TLine* fcllLine = new TLine(fcll, 0, fcll, 1); TLine* fculLine = new TLine(fcul, 0, fcul, 1); fcllLine->SetLineColor(kRed); fculLine->SetLineColor(kRed); fcllLine->Draw("same"); fculLine->Draw("same"); dataCanvas->Update(); } // Plot MCMC interval and print some statistics MCMCIntervalPlot mcPlot(*mcInt); mcPlot.SetLineColor(kMagenta); mcPlot.SetLineWidth(2); mcPlot.Draw("same"); double mcul = mcInt->UpperLimit(*s); double mcll = mcInt->LowerLimit(*s); cout << "MCMC lower limit on s = " << mcll << endl; cout << "MCMC upper limit on s = " << mcul << endl; cout << "MCMC Actual confidence level: " << mcInt->GetActualConfidenceLevel() << endl; // 3-d plot of the parameter points dataCanvas->cd(2); // also plot the points in the markov chain TTree& chain = ((RooTreeDataStore*) mcInt->GetChainAsDataSet()->store())->tree(); chain.SetMarkerStyle(6); chain.SetMarkerColor(kRed); chain.Draw("s:ratioSigEff:ratioBkgEff","weight_MarkovChain_local_","box"); // 3-d box proporional to posterior // the points used in the profile construction TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree(); parameterScan.SetMarkerStyle(24); parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same"); delete wspace; delete lrint; delete mcInt; delete fcint; delete data; /// print timing info t.Stop(); t.Print(); } rs101_limitexample.C:1 rs101_limitexample.C:2 rs101_limitexample.C:3 rs101_limitexample.C:4 rs101_limitexample.C:5 rs101_limitexample.C:6 rs101_limitexample.C:7 rs101_limitexample.C:8 rs101_limitexample.C:9 rs101_limitexample.C:10 rs101_limitexample.C:11 rs101_limitexample.C:12 rs101_limitexample.C:13 rs101_limitexample.C:14 rs101_limitexample.C:15 rs101_limitexample.C:16 rs101_limitexample.C:17 rs101_limitexample.C:18 rs101_limitexample.C:19 rs101_limitexample.C:20 rs101_limitexample.C:21 rs101_limitexample.C:22 rs101_limitexample.C:23 rs101_limitexample.C:24 rs101_limitexample.C:25 rs101_limitexample.C:26 rs101_limitexample.C:27 rs101_limitexample.C:28 rs101_limitexample.C:29 rs101_limitexample.C:30 rs101_limitexample.C:31 rs101_limitexample.C:32 rs101_limitexample.C:33 rs101_limitexample.C:34 rs101_limitexample.C:35 rs101_limitexample.C:36 rs101_limitexample.C:37 rs101_limitexample.C:38 rs101_limitexample.C:39 rs101_limitexample.C:40 rs101_limitexample.C:41 rs101_limitexample.C:42 rs101_limitexample.C:43 rs101_limitexample.C:44 rs101_limitexample.C:45 rs101_limitexample.C:46 rs101_limitexample.C:47 rs101_limitexample.C:48 rs101_limitexample.C:49 rs101_limitexample.C:50 rs101_limitexample.C:51 rs101_limitexample.C:52 rs101_limitexample.C:53 rs101_limitexample.C:54 rs101_limitexample.C:55 rs101_limitexample.C:56 rs101_limitexample.C:57 rs101_limitexample.C:58 rs101_limitexample.C:59 rs101_limitexample.C:60 rs101_limitexample.C:61 rs101_limitexample.C:62 rs101_limitexample.C:63 rs101_limitexample.C:64 rs101_limitexample.C:65 rs101_limitexample.C:66 rs101_limitexample.C:67 rs101_limitexample.C:68 rs101_limitexample.C:69 rs101_limitexample.C:70 rs101_limitexample.C:71 rs101_limitexample.C:72 rs101_limitexample.C:73 rs101_limitexample.C:74 rs101_limitexample.C:75 rs101_limitexample.C:76 rs101_limitexample.C:77 rs101_limitexample.C:78 rs101_limitexample.C:79 rs101_limitexample.C:80 rs101_limitexample.C:81 rs101_limitexample.C:82 rs101_limitexample.C:83 rs101_limitexample.C:84 rs101_limitexample.C:85 rs101_limitexample.C:86 rs101_limitexample.C:87 rs101_limitexample.C:88 rs101_limitexample.C:89 rs101_limitexample.C:90 rs101_limitexample.C:91 rs101_limitexample.C:92 rs101_limitexample.C:93 rs101_limitexample.C:94 rs101_limitexample.C:95 rs101_limitexample.C:96 rs101_limitexample.C:97 rs101_limitexample.C:98 rs101_limitexample.C:99 rs101_limitexample.C:100 rs101_limitexample.C:101 rs101_limitexample.C:102 rs101_limitexample.C:103 rs101_limitexample.C:104 rs101_limitexample.C:105 rs101_limitexample.C:106 rs101_limitexample.C:107 rs101_limitexample.C:108 rs101_limitexample.C:109 rs101_limitexample.C:110 rs101_limitexample.C:111 rs101_limitexample.C:112 rs101_limitexample.C:113 rs101_limitexample.C:114 rs101_limitexample.C:115 rs101_limitexample.C:116 rs101_limitexample.C:117 rs101_limitexample.C:118 rs101_limitexample.C:119 rs101_limitexample.C:120 rs101_limitexample.C:121 rs101_limitexample.C:122 rs101_limitexample.C:123 rs101_limitexample.C:124 rs101_limitexample.C:125 rs101_limitexample.C:126 rs101_limitexample.C:127 rs101_limitexample.C:128 rs101_limitexample.C:129 rs101_limitexample.C:130 rs101_limitexample.C:131 rs101_limitexample.C:132 rs101_limitexample.C:133 rs101_limitexample.C:134 rs101_limitexample.C:135 rs101_limitexample.C:136 rs101_limitexample.C:137 rs101_limitexample.C:138 rs101_limitexample.C:139 rs101_limitexample.C:140 rs101_limitexample.C:141 rs101_limitexample.C:142 rs101_limitexample.C:143 rs101_limitexample.C:144 rs101_limitexample.C:145 rs101_limitexample.C:146 rs101_limitexample.C:147 rs101_limitexample.C:148 rs101_limitexample.C:149 rs101_limitexample.C:150 rs101_limitexample.C:151 rs101_limitexample.C:152 rs101_limitexample.C:153 rs101_limitexample.C:154 rs101_limitexample.C:155 rs101_limitexample.C:156 rs101_limitexample.C:157 rs101_limitexample.C:158 rs101_limitexample.C:159 rs101_limitexample.C:160 rs101_limitexample.C:161 rs101_limitexample.C:162 rs101_limitexample.C:163 rs101_limitexample.C:164 rs101_limitexample.C:165 rs101_limitexample.C:166 rs101_limitexample.C:167 rs101_limitexample.C:168 rs101_limitexample.C:169 rs101_limitexample.C:170 rs101_limitexample.C:171 rs101_limitexample.C:172 rs101_limitexample.C:173 rs101_limitexample.C:174 rs101_limitexample.C:175 rs101_limitexample.C:176 rs101_limitexample.C:177 rs101_limitexample.C:178 rs101_limitexample.C:179 rs101_limitexample.C:180 rs101_limitexample.C:181 rs101_limitexample.C:182 rs101_limitexample.C:183 rs101_limitexample.C:184 rs101_limitexample.C:185 rs101_limitexample.C:186 rs101_limitexample.C:187 rs101_limitexample.C:188 rs101_limitexample.C:189 rs101_limitexample.C:190 rs101_limitexample.C:191 rs101_limitexample.C:192 rs101_limitexample.C:193 rs101_limitexample.C:194 rs101_limitexample.C:195 rs101_limitexample.C:196 rs101_limitexample.C:197 rs101_limitexample.C:198 rs101_limitexample.C:199 rs101_limitexample.C:200 rs101_limitexample.C:201 rs101_limitexample.C:202 rs101_limitexample.C:203 rs101_limitexample.C:204 |
|