rs301_splot.C: SPlot tutorial | Roostats tutorials | rs401d_FeldmanCousins.C: 'Neutrino Oscillation Example from Feldman & Cousins' |
///////////////////////////////////////////////////////////////////////// // // Produces an interval on the mean signal in a number counting // experiment with known background using the Feldman-Cousins technique. // date Jan. 2009 // updated June 2010 // // Using the RooStats FeldmanCousins tool with 200 bins // it takes 1 min and the interval is [0.2625, 10.6125] // with a step size of 0.075. // The interval in Feldman & Cousins's original paper is [.29, 10.81] // Phys.Rev.D57:3873-3889,1998. ///////////////////////////////////////////////////////////////////////// #include "RooGlobalFunc.h" #include "RooStats/ConfInterval.h" #include "RooStats/PointSetInterval.h" #include "RooStats/ConfidenceBelt.h" #include "RooStats/FeldmanCousins.h" #include "RooStats/ModelConfig.h" #include "RooWorkspace.h" #include "RooDataSet.h" #include "RooRealVar.h" #include "RooConstVar.h" #include "RooAddition.h" #include "RooDataHist.h" #include "RooPoisson.h" #include "RooPlot.h" #include "TCanvas.h" #include "TTree.h" #include "TH1F.h" #include "TMarker.h" #include "TStopwatch.h" #include <iostream> // use this order for safety on library loading using namespace RooFit ; using namespace RooStats ; void rs401c_FeldmanCousins() { // to time the macro... about 30 s TStopwatch t; t.Start(); // make a simple model RooRealVar x("x","", 1,0,50); RooRealVar mu("mu","", 2.5,0, 15); // with a limit on mu>=0 RooConstVar b("b","", 3.); RooAddition mean("mean","",RooArgList(mu,b)); RooPoisson pois("pois", "", x, mean); RooArgSet parameters(mu); // create a toy dataset RooDataSet* data = pois.generate(RooArgSet(x), 1); data->Print("v"); TCanvas* dataCanvas = new TCanvas("dataCanvas"); RooPlot* frame = x.frame(); data->plotOn(frame); frame->Draw(); dataCanvas->Update(); RooWorkspace* w = new RooWorkspace(); ModelConfig modelConfig("poissonProblem",w); modelConfig.SetPdf(pois); modelConfig.SetParametersOfInterest(parameters); modelConfig.SetObservables(RooArgSet(x)); w->Print(); //////// show use of Feldman-Cousins RooStats::FeldmanCousins fc(*data,modelConfig); fc.SetTestSize(.05); // set size of test 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 // use the Feldman-Cousins tool PointSetInterval* interval = (PointSetInterval*)fc.GetInterval(); // make a canvas for plots TCanvas* intervalCanvas = new TCanvas("intervalCanvas"); std::cout << "is this point in the interval? " << interval->IsInInterval(parameters) << std::endl; std::cout << "interval is ["<< interval->LowerLimit(mu) << ", " << interval->UpperLimit(mu) << "]" << endl; // using 200 bins it takes 1 min and the answer is // interval is [0.2625, 10.6125] with a step size of .075 // The interval in Feldman & Cousins's original paper is [.29, 10.81] // Phys.Rev.D57:3873-3889,1998. // No dedicated plotting class yet, so do it by hand: RooDataHist* parameterScan = (RooDataHist*) fc.GetPointsToScan(); TH1F* hist = (TH1F*) parameterScan->createHistogram("mu",30); hist->Draw(); RooArgSet* tmpPoint; // loop over points to test for(Int_t i=0; i<parameterScan->numEntries(); ++i){ // cout << "on parameter point " << i << " out of " << parameterScan->numEntries() << endl; // get a parameter point from the list of points to test. tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp"); TMarker* mark = new TMarker(tmpPoint->getRealValue("mu"), 1, 25); if (interval->IsInInterval( *tmpPoint ) ) mark->SetMarkerColor(kBlue); else mark->SetMarkerColor(kRed); mark->Draw("s"); //delete tmpPoint; // delete mark; } t.Stop(); t.Print(); } rs401c_FeldmanCousins.C:1 rs401c_FeldmanCousins.C:2 rs401c_FeldmanCousins.C:3 rs401c_FeldmanCousins.C:4 rs401c_FeldmanCousins.C:5 rs401c_FeldmanCousins.C:6 rs401c_FeldmanCousins.C:7 rs401c_FeldmanCousins.C:8 rs401c_FeldmanCousins.C:9 rs401c_FeldmanCousins.C:10 rs401c_FeldmanCousins.C:11 rs401c_FeldmanCousins.C:12 rs401c_FeldmanCousins.C:13 rs401c_FeldmanCousins.C:14 rs401c_FeldmanCousins.C:15 rs401c_FeldmanCousins.C:16 rs401c_FeldmanCousins.C:17 rs401c_FeldmanCousins.C:18 rs401c_FeldmanCousins.C:19 rs401c_FeldmanCousins.C:20 rs401c_FeldmanCousins.C:21 rs401c_FeldmanCousins.C:22 rs401c_FeldmanCousins.C:23 rs401c_FeldmanCousins.C:24 rs401c_FeldmanCousins.C:25 rs401c_FeldmanCousins.C:26 rs401c_FeldmanCousins.C:27 rs401c_FeldmanCousins.C:28 rs401c_FeldmanCousins.C:29 rs401c_FeldmanCousins.C:30 rs401c_FeldmanCousins.C:31 rs401c_FeldmanCousins.C:32 rs401c_FeldmanCousins.C:33 rs401c_FeldmanCousins.C:34 rs401c_FeldmanCousins.C:35 rs401c_FeldmanCousins.C:36 rs401c_FeldmanCousins.C:37 rs401c_FeldmanCousins.C:38 rs401c_FeldmanCousins.C:39 rs401c_FeldmanCousins.C:40 rs401c_FeldmanCousins.C:41 rs401c_FeldmanCousins.C:42 rs401c_FeldmanCousins.C:43 rs401c_FeldmanCousins.C:44 rs401c_FeldmanCousins.C:45 rs401c_FeldmanCousins.C:46 rs401c_FeldmanCousins.C:47 rs401c_FeldmanCousins.C:48 rs401c_FeldmanCousins.C:49 rs401c_FeldmanCousins.C:50 rs401c_FeldmanCousins.C:51 rs401c_FeldmanCousins.C:52 rs401c_FeldmanCousins.C:53 rs401c_FeldmanCousins.C:54 rs401c_FeldmanCousins.C:55 rs401c_FeldmanCousins.C:56 rs401c_FeldmanCousins.C:57 rs401c_FeldmanCousins.C:58 rs401c_FeldmanCousins.C:59 rs401c_FeldmanCousins.C:60 rs401c_FeldmanCousins.C:61 rs401c_FeldmanCousins.C:62 rs401c_FeldmanCousins.C:63 rs401c_FeldmanCousins.C:64 rs401c_FeldmanCousins.C:65 rs401c_FeldmanCousins.C:66 rs401c_FeldmanCousins.C:67 rs401c_FeldmanCousins.C:68 rs401c_FeldmanCousins.C:69 rs401c_FeldmanCousins.C:70 rs401c_FeldmanCousins.C:71 rs401c_FeldmanCousins.C:72 rs401c_FeldmanCousins.C:73 rs401c_FeldmanCousins.C:74 rs401c_FeldmanCousins.C:75 rs401c_FeldmanCousins.C:76 rs401c_FeldmanCousins.C:77 rs401c_FeldmanCousins.C:78 rs401c_FeldmanCousins.C:79 rs401c_FeldmanCousins.C:80 rs401c_FeldmanCousins.C:81 rs401c_FeldmanCousins.C:82 rs401c_FeldmanCousins.C:83 rs401c_FeldmanCousins.C:84 rs401c_FeldmanCousins.C:85 rs401c_FeldmanCousins.C:86 rs401c_FeldmanCousins.C:87 rs401c_FeldmanCousins.C:88 rs401c_FeldmanCousins.C:89 rs401c_FeldmanCousins.C:90 rs401c_FeldmanCousins.C:91 rs401c_FeldmanCousins.C:92 rs401c_FeldmanCousins.C:93 rs401c_FeldmanCousins.C:94 rs401c_FeldmanCousins.C:95 rs401c_FeldmanCousins.C:96 rs401c_FeldmanCousins.C:97 rs401c_FeldmanCousins.C:98 rs401c_FeldmanCousins.C:99 rs401c_FeldmanCousins.C:100 rs401c_FeldmanCousins.C:101 rs401c_FeldmanCousins.C:102 rs401c_FeldmanCousins.C:103 rs401c_FeldmanCousins.C:104 rs401c_FeldmanCousins.C:105 rs401c_FeldmanCousins.C:106 rs401c_FeldmanCousins.C:107 rs401c_FeldmanCousins.C:108 rs401c_FeldmanCousins.C:109 rs401c_FeldmanCousins.C:110 rs401c_FeldmanCousins.C:111 rs401c_FeldmanCousins.C:112 rs401c_FeldmanCousins.C:113 rs401c_FeldmanCousins.C:114 rs401c_FeldmanCousins.C:115 rs401c_FeldmanCousins.C:116 rs401c_FeldmanCousins.C:117 rs401c_FeldmanCousins.C:118 rs401c_FeldmanCousins.C:119 rs401c_FeldmanCousins.C:120 rs401c_FeldmanCousins.C:121 rs401c_FeldmanCousins.C:122 rs401c_FeldmanCousins.C:123 rs401c_FeldmanCousins.C:124 rs401c_FeldmanCousins.C:125 rs401c_FeldmanCousins.C:126 rs401c_FeldmanCousins.C:127 rs401c_FeldmanCousins.C:128 rs401c_FeldmanCousins.C:129 rs401c_FeldmanCousins.C:130 rs401c_FeldmanCousins.C:131 |
|