ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
OneSidedFrequentistUpperLimitWithBands.C
Go to the documentation of this file.
1 /*
2 OneSidedFrequentistUpperLimitWithBands
3 
4 Author: Kyle Cranmer,
5 Contributions from Haichen Wang and Daniel Whiteson
6 date: Dec. 2010 - Feb. 2011
7 v1. Jan 28
8 
9 This is a standard demo that can be used with any ROOT file
10 prepared in the standard way. You specify:
11  - name for input ROOT file
12  - name of workspace inside ROOT file that holds model and data
13  - name of ModelConfig that specifies details for calculator tools
14  - name of dataset
15 
16 With default parameters the macro will attempt to run the
17 standard hist2workspace example and read the ROOT file
18 that it produces.
19 
20 The first ~100 lines define a new test statistic, then the main macro starts.
21 You may want to control:
22  double confidenceLevel=0.95;
23  int nPointsToScan = 30;
24  int nToyMC = 200;
25 
26 This uses a modified version of the profile likelihood ratio as
27 a test statistic for upper limits (eg. test stat = 0 if muhat>mu).
28 
29 Based on the observed data, one defines a set of parameter points
30 to be tested based on the value of the parameter of interest
31 and the conditional MLE (eg. profiled) values of the nuisance parameters.
32 
33 At each parameter point, pseudo-experiments are generated using this
34 fixed reference model and then the test statistic is evaluated.
35 Note, the nuisance parameters are floating in the fits. For each point,
36 the threshold that defines the 95% acceptance region is found. This
37 forms a "Confidence Belt".
38 
39 After constructing the confidence belt, one can find the confidence
40 interval for any particular dataset by finding the intersection
41 of the observed test statistic and the confidence belt. First
42 this is done on the observed data to get an observed 1-sided upper limt.
43 
44 Finally, there expected limit and bands (from background-only) are
45 formed by generating background-only data and finding the upper limit.
46 This is done by hand for now, will later be part of the RooStats tools.
47 
48 On a technical note, this technique is NOT the Feldman-Cousins technique,
49 because that is a 2-sided interval BY DEFINITION. However, like the
50 Feldman-Cousins technique this is a Neyman-Construction. For technical
51 reasons the easiest way to implement this right now is to use the
52 FeldmanCousins tool and then change the test statistic that it is using.
53 
54 Building the confidence belt can be computationally expensive. Once it is built,
55 one could save it to a file and use it in a separate step.
56 
57 We can use PROOF to speed things along in parallel, however,
58 the test statistic has to be installed on the workers
59 so either turn off PROOF or include the modified test statistic
60 in your $ROOTSYS/roofit/roostats/inc directory,
61 add the additional line to the LinkDef.h file,
62 and recompile root.
63 
64 Note, if you have a boundary on the parameter of interest (eg. cross-section)
65 the threshold on the one-sided test statistic starts off very small because we
66 are only including downward fluctuations. You can see the threshold in these printouts:
67 
68 [#0] PROGRESS:Generation -- generated toys: 500 / 999
69 NeymanConstruction: Prog: 12/50 total MC = 39 this test stat = 0
70  SigXsecOverSM=0.69 alpha_syst1=0.136515 alpha_syst3=0.425415 beta_syst2=1.08496 [-1e+30, 0.011215] in interval = 1
71 
72 this tells you the values of the parameters being used to generate the pseudo-experiments
73 and the threshold in this case is 0.011215. One would expect for 95% that the threshold
74 would be ~1.35 once the cross-section is far enough away from 0 that it is essentially
75 unaffected by the boundary. As one reaches the last points in the scan, the
76 theshold starts to get artificially high. This is because the range of the parameter in
77 the fit is the same as the range in the scan. In the future, these should be independently
78 controled, but they are not now. As a result the ~50% of pseudo-experiments that have an
79 upward fluctuation end up with muhat = muMax. Because of this, the upper range of the
80 parameter should be well above the expected upper limit... but not too high or one will
81 need a very large value of nPointsToScan to resolve the relevant region. This can be
82 improved, but this is the first version of this script.
83 
84 Important note: when the model includes external constraint terms, like a Gaussian
85 constraint to a nuisance parameter centered around some nominal value there is
86 a subtlety. The asymptotic results are all based on the assumption that all the
87 measurements fluctuate... including the nominal values from auxiliary measurements.
88 If these do not fluctuate, this corresponds to an "conditional ensemble". The
89 result is that the distribution of the test statistic can become very non-chi^2.
90 This results in thresholds that become very large. This can be seen in the following
91 thought experiment. Say the model is
92  Pois(N | s + b) * G(b0|b,sigma)
93 where G(b0|b,sigma) is the external constraint and b0 is 100. If N is also 100
94 then the profiled value of b given s is going to be some trade off betwen 100-s and b0.
95 If sigma is \sqrt(N), then the profiled value of b is probably 100 - s/2 So for
96 s=60 we are going to have a profiled value of b~70. Now when we generate pseudo-experiments
97 for s=60, b=70 we will have N~130 and the average shat will be 30, not 60. In practice,
98 this is only an issue for values of s that are very excluded. For values of s near the 95%
99 limit this should not be a big effect. This can be avoided if the nominal values of the constraints also fluctuate, but that requires that those parameters are RooRealVars in the model.
100 This version does not deal with this issue, but it will be addressed in a future version.
101 */
102 
103 #include "TFile.h"
104 #include "TROOT.h"
105 #include "TH1F.h"
106 #include "TCanvas.h"
107 #include "TSystem.h"
108 
109 #include "RooWorkspace.h"
110 #include "RooSimultaneous.h"
111 #include "RooAbsData.h"
112 
113 #include "RooStats/ModelConfig.h"
114 #include "RooStats/FeldmanCousins.h"
115 #include "RooStats/ToyMCSampler.h"
117 #include "RooStats/ConfidenceBelt.h"
118 
119 #include "RooStats/RooStatsUtils.h"
121 
122 using namespace RooFit;
123 using namespace RooStats;
124 
125 bool useProof = false; // flag to control whether to use Proof
126 int nworkers = 0; // number of workers (default use all available cores)
127 
128 /////////////////////////////////////////////////////////////////////////
129 // The actual macro
130 
132  const char* workspaceName = "combined",
133  const char* modelConfigName = "ModelConfig",
134  const char* dataName = "obsData") {
135 
136 
137 #ifdef __CINT__
138  cout << "DO NOT RUN WITH CINT: we are using a custom test statistic ";
139  cout << "which requires that this tutorial must be compiled ";
140  cout << "with ACLIC" << endl;
141  return;
142 #endif
143 
144  double confidenceLevel=0.95;
145  int nPointsToScan = 20;
146  int nToyMC = 200;
147 
148  /////////////////////////////////////////////////////////////
149  // First part is just to access a user-defined file
150  // or create the standard example file if it doesn't exist
151  ////////////////////////////////////////////////////////////
152  const char* filename = "";
153  if (!strcmp(infile,"")) {
154  filename = "results/example_combined_GaussExample_model.root";
155  bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
156  // if file does not exists generate with histfactory
157  if (!fileExist) {
158 #ifdef _WIN32
159  cout << "HistFactory file cannot be generated on Windows - exit" << endl;
160  return;
161 #endif
162  // Normally this would be run on the command line
163  cout <<"will run standard hist2workspace example"<<endl;
164  gROOT->ProcessLine(".! prepareHistFactory .");
165  gROOT->ProcessLine(".! hist2workspace config/example.xml");
166  cout <<"\n\n---------------------"<<endl;
167  cout <<"Done creating example input"<<endl;
168  cout <<"---------------------\n\n"<<endl;
169  }
170 
171  }
172  else
173  filename = infile;
174 
175  // Try to open the file
176  TFile *file = TFile::Open(filename);
177 
178  // if input file was specified byt not found, quit
179  if(!file ){
180  cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
181  return;
182  }
183 
184 
185  /////////////////////////////////////////////////////////////
186  // Now get the data and workspace
187  ////////////////////////////////////////////////////////////
188 
189  // get the workspace out of the file
190  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
191  if(!w){
192  cout <<"workspace not found" << endl;
193  return;
194  }
195 
196  // get the modelConfig out of the file
197  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);
198 
199  // get the modelConfig out of the file
200  RooAbsData* data = w->data(dataName);
201 
202  // make sure ingredients are found
203  if(!data || !mc){
204  w->Print();
205  cout << "data or ModelConfig was not found" <<endl;
206  return;
207  }
208 
209  /////////////////////////////////////////////////////////////
210  // Now get the POI for convenience
211  // you may want to adjust the range of your POI
212  ////////////////////////////////////////////////////////////
213  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
214  // firstPOI->setMin(0);
215  // firstPOI->setMax(10);
216 
217  /////////////////////////////////////////////
218  // create and use the FeldmanCousins tool
219  // to find and plot the 95% confidence interval
220  // on the parameter of interest as specified
221  // in the model config
222  // REMEMBER, we will change the test statistic
223  // so this is NOT a Feldman-Cousins interval
224  FeldmanCousins fc(*data,*mc);
225  fc.SetConfidenceLevel(confidenceLevel);
226  // fc.AdditionalNToysFactor(0.25); // degrade/improve sampling that defines confidence belt
227  // fc.UseAdaptiveSampling(true); // speed it up a bit, don't use for expectd limits
228  fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
229  fc.CreateConfBelt(true); // save the information in the belt for plotting
230 
231  /////////////////////////////////////////////
232  // Feldman-Cousins is a unified limit by definition
233  // but the tool takes care of a few things for us like which values
234  // of the nuisance parameters should be used to generate toys.
235  // so let's just change the test statistic and realize this is
236  // no longer "Feldman-Cousins" but is a fully frequentist Neyman-Construction.
237  // ProfileLikelihoodTestStatModified onesided(*mc->GetPdf());
238  // fc.GetTestStatSampler()->SetTestStatistic(&onesided);
239  // ((ToyMCSampler*) fc.GetTestStatSampler())->SetGenerateBinned(true);
240  ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
241  ProfileLikelihoodTestStat* testStat = dynamic_cast<ProfileLikelihoodTestStat*>(toymcsampler->GetTestStatistic());
242  testStat->SetOneSided(true);
243 
244  // Since this tool needs to throw toy MC the PDF needs to be
245  // extended or the tool needs to know how many entries in a dataset
246  // per pseudo experiment.
247  // In the 'number counting form' where the entries in the dataset
248  // are counts, and not values of discriminating variables, the
249  // datasets typically only have one entry and the PDF is not
250  // extended.
251  if(!mc->GetPdf()->canBeExtended()){
252  if(data->numEntries()==1)
253  fc.FluctuateNumDataEntries(false);
254  else
255  cout <<"Not sure what to do about this model" <<endl;
256  }
257 
258  // We can use PROOF to speed things along in parallel
259  // However, the test statistic has to be installed on the workers
260  // so either turn off PROOF or include the modified test statistic
261  // in your $ROOTSYS/roofit/roostats/inc directory,
262  // add the additional line to the LinkDef.h file,
263  // and recompile root.
264  if (useProof) {
265  ProofConfig pc(*w, nworkers, "", false);
266  toymcsampler->SetProofConfig(&pc); // enable proof
267  }
268 
269  if(mc->GetGlobalObservables()){
270  cout << "will use global observables for unconditional ensemble"<<endl;
271  mc->GetGlobalObservables()->Print();
272  toymcsampler->SetGlobalObservables(*mc->GetGlobalObservables());
273  }
274 
275 
276  // Now get the interval
277  PointSetInterval* interval = fc.GetInterval();
278  ConfidenceBelt* belt = fc.GetConfidenceBelt();
279 
280  // print out the iterval on the first Parameter of Interest
281  cout << "\n95% interval on " <<firstPOI->GetName()<<" is : ["<<
282  interval->LowerLimit(*firstPOI) << ", "<<
283  interval->UpperLimit(*firstPOI) <<"] "<<endl;
284 
285  // get observed UL and value of test statistic evaluated there
286  RooArgSet tmpPOI(*firstPOI);
287  double observedUL = interval->UpperLimit(*firstPOI);
288  firstPOI->setVal(observedUL);
289  double obsTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*data,tmpPOI);
290 
291 
292  // Ask the calculator which points were scanned
293  RooDataSet* parameterScan = (RooDataSet*) fc.GetPointsToScan();
294  RooArgSet* tmpPoint;
295 
296  // make a histogram of parameter vs. threshold
297  TH1F* histOfThresholds = new TH1F("histOfThresholds","",
298  parameterScan->numEntries(),
299  firstPOI->getMin(),
300  firstPOI->getMax());
301  histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
302  histOfThresholds->GetYaxis()->SetTitle("Threshold");
303 
304  // loop through the points that were tested and ask confidence belt
305  // what the upper/lower thresholds were.
306  // For FeldmanCousins, the lower cut off is always 0
307  for(Int_t i=0; i<parameterScan->numEntries(); ++i){
308  tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
309  //cout <<"get threshold"<<endl;
310  double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
311  double poiVal = tmpPoint->getRealValue(firstPOI->GetName()) ;
312  histOfThresholds->Fill(poiVal,arMax);
313  }
314  TCanvas* c1 = new TCanvas();
315  c1->Divide(2);
316  c1->cd(1);
317  histOfThresholds->SetMinimum(0);
318  histOfThresholds->Draw();
319  c1->cd(2);
320 
321  /////////////////////////////////////////////////////////////
322  // Now we generate the expected bands and power-constriant
323  ////////////////////////////////////////////////////////////
324 
325  // First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
326  RooAbsReal* nll = mc->GetPdf()->createNLL(*data);
327  RooAbsReal* profile = nll->createProfile(*mc->GetParametersOfInterest());
328  firstPOI->setVal(0.);
329  profile->getVal(); // this will do fit and set nuisance parameters to profiled values
330  RooArgSet* poiAndNuisance = new RooArgSet();
331  if(mc->GetNuisanceParameters())
332  poiAndNuisance->add(*mc->GetNuisanceParameters());
333  poiAndNuisance->add(*mc->GetParametersOfInterest());
334  w->saveSnapshot("paramsToGenerateData",*poiAndNuisance);
335  RooArgSet* paramsToGenerateData = (RooArgSet*) poiAndNuisance->snapshot();
336  cout << "\nWill use these parameter points to generate pseudo data for bkg only" << endl;
337  paramsToGenerateData->Print("v");
338 
339 
340  RooArgSet unconditionalObs;
341  unconditionalObs.add(*mc->GetObservables());
342  unconditionalObs.add(*mc->GetGlobalObservables()); // comment this out for the original conditional ensemble
343 
344  double CLb=0;
345  double CLbinclusive=0;
346 
347  // Now we generate background only and find distribution of upper limits
348  TH1F* histOfUL = new TH1F("histOfUL","",100,0,firstPOI->getMax());
349  histOfUL->GetXaxis()->SetTitle("Upper Limit (background only)");
350  histOfUL->GetYaxis()->SetTitle("Entries");
351  for(int imc=0; imc<nToyMC; ++imc){
352 
353  // set parameters back to values for generating pseudo data
354  // cout << "\n get current nuis, set vals, print again" << endl;
355  w->loadSnapshot("paramsToGenerateData");
356  // poiAndNuisance->Print("v");
357 
358  RooDataSet* toyData = 0;
359  // now generate a toy dataset
360  if(!mc->GetPdf()->canBeExtended()){
361  if(data->numEntries()==1)
362  toyData = mc->GetPdf()->generate(*mc->GetObservables(),1);
363  else
364  cout <<"Not sure what to do about this model" <<endl;
365  } else{
366  // cout << "generating extended dataset"<<endl;
367  toyData = mc->GetPdf()->generate(*mc->GetObservables(),Extended());
368  }
369 
370  // generate global observables
371  // need to be careful for simpdf
372  // RooDataSet* globalData = mc->GetPdf()->generate(*mc->GetGlobalObservables(),1);
373 
374  RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>(mc->GetPdf());
375  if(!simPdf){
376  RooDataSet *one = mc->GetPdf()->generate(*mc->GetGlobalObservables(), 1);
377  const RooArgSet *values = one->get();
378  RooArgSet *allVars = mc->GetPdf()->getVariables();
379  *allVars = *values;
380  delete allVars;
381  delete values;
382  delete one;
383  } else {
384 
385  //try fix for sim pdf
386  TIterator* iter = simPdf->indexCat().typeIterator() ;
387  RooCatType* tt = NULL;
388  while((tt=(RooCatType*) iter->Next())) {
389 
390  // Get pdf associated with state from simpdf
391  RooAbsPdf* pdftmp = simPdf->getPdf(tt->GetName()) ;
392 
393  // Generate only global variables defined by the pdf associated with this state
394  RooArgSet* globtmp = pdftmp->getObservables(*mc->GetGlobalObservables()) ;
395  RooDataSet* tmp = pdftmp->generate(*globtmp,1) ;
396 
397  // Transfer values to output placeholder
398  *globtmp = *tmp->get(0) ;
399 
400  // Cleanup
401  delete globtmp ;
402  delete tmp ;
403  }
404  }
405 
406  // globalData->Print("v");
407  // unconditionalObs = *globalData->get();
408  // mc->GetGlobalObservables()->Print("v");
409  // delete globalData;
410  // cout << "toy data = " << endl;
411  // toyData->get()->Print("v");
412 
413  // get test stat at observed UL in observed data
414  firstPOI->setVal(observedUL);
415  double toyTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);
416  // toyData->get()->Print("v");
417  // cout <<"obsTSatObsUL " <<obsTSatObsUL << "toyTS " << toyTSatObsUL << endl;
418  if(obsTSatObsUL < toyTSatObsUL) // not sure about <= part yet
419  CLb+= (1.)/nToyMC;
420  if(obsTSatObsUL <= toyTSatObsUL) // not sure about <= part yet
421  CLbinclusive+= (1.)/nToyMC;
422 
423 
424  // loop over points in belt to find upper limit for this toy data
425  double thisUL = 0;
426  for(Int_t i=0; i<parameterScan->numEntries(); ++i){
427  tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
428  double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
429  firstPOI->setVal( tmpPoint->getRealValue(firstPOI->GetName()) );
430  // double thisTS = profile->getVal();
431  double thisTS = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);
432 
433  // cout << "poi = " << firstPOI->getVal()
434  // << " max is " << arMax << " this profile = " << thisTS << endl;
435  // cout << "thisTS = " << thisTS<<endl;
436  if(thisTS<=arMax){
437  thisUL = firstPOI->getVal();
438  } else{
439  break;
440  }
441  }
442 
443 
444 
445  /*
446  // loop over points in belt to find upper limit for this toy data
447  double thisUL = 0;
448  for(Int_t i=0; i<histOfThresholds->GetNbinsX(); ++i){
449  tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
450  cout <<"---------------- "<<i<<endl;
451  tmpPoint->Print("v");
452  cout << "from hist " << histOfThresholds->GetBinCenter(i+1) <<endl;
453  double arMax = histOfThresholds->GetBinContent(i+1);
454  // cout << " threhold from Hist = aMax " << arMax<<endl;
455  // double arMax2 = belt->GetAcceptanceRegionMax(*tmpPoint);
456  // cout << "from scan arMax2 = "<< arMax2 << endl; // not the same due to TH1F not TH1D
457  // cout << "scan - hist" << arMax2-arMax << endl;
458  firstPOI->setVal( histOfThresholds->GetBinCenter(i+1));
459  // double thisTS = profile->getVal();
460  double thisTS = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);
461 
462  // cout << "poi = " << firstPOI->getVal()
463  // << " max is " << arMax << " this profile = " << thisTS << endl;
464  // cout << "thisTS = " << thisTS<<endl;
465 
466  // NOTE: need to add a small epsilon term for single precision vs. double precision
467  if(thisTS<=arMax + 1e-7){
468  thisUL = firstPOI->getVal();
469  } else{
470  break;
471  }
472  }
473  */
474 
475  histOfUL->Fill(thisUL);
476 
477  // for few events, data is often the same, and UL is often the same
478  // cout << "thisUL = " << thisUL<<endl;
479 
480  delete toyData;
481  }
482  histOfUL->Draw();
483  c1->SaveAs("one-sided_upper_limit_output.pdf");
484 
485  // if you want to see a plot of the sampling distribution for a particular scan point:
486  /*
487  SamplingDistPlot sampPlot;
488  int indexInScan = 0;
489  tmpPoint = (RooArgSet*) parameterScan->get(indexInScan)->clone("temp");
490  firstPOI->setVal( tmpPoint->getRealValue(firstPOI->GetName()) );
491  toymcsampler->SetParametersForTestStat(tmpPOI);
492  SamplingDistribution* samp = toymcsampler->GetSamplingDistribution(*tmpPoint);
493  sampPlot.AddSamplingDistribution(samp);
494  sampPlot.Draw();
495  */
496 
497  // Now find bands and power constraint
498  Double_t* bins = histOfUL->GetIntegral();
499  TH1F* cumulative = (TH1F*) histOfUL->Clone("cumulative");
500  cumulative->SetContent(bins);
501  double band2sigDown, band1sigDown, bandMedian, band1sigUp,band2sigUp;
502  for(int i=1; i<=cumulative->GetNbinsX(); ++i){
503  if(bins[i]<RooStats::SignificanceToPValue(2))
504  band2sigDown=cumulative->GetBinCenter(i);
505  if(bins[i]<RooStats::SignificanceToPValue(1))
506  band1sigDown=cumulative->GetBinCenter(i);
507  if(bins[i]<0.5)
508  bandMedian=cumulative->GetBinCenter(i);
509  if(bins[i]<RooStats::SignificanceToPValue(-1))
510  band1sigUp=cumulative->GetBinCenter(i);
511  if(bins[i]<RooStats::SignificanceToPValue(-2))
512  band2sigUp=cumulative->GetBinCenter(i);
513  }
514  cout << "-2 sigma band " << band2sigDown << endl;
515  cout << "-1 sigma band " << band1sigDown << " [Power Constriant)]" << endl;
516  cout << "median of band " << bandMedian << endl;
517  cout << "+1 sigma band " << band1sigUp << endl;
518  cout << "+2 sigma band " << band2sigUp << endl;
519 
520  // print out the iterval on the first Parameter of Interest
521  cout << "\nobserved 95% upper-limit "<< interval->UpperLimit(*firstPOI) <<endl;
522  cout << "CLb strict [P(toy>obs|0)] for observed 95% upper-limit "<< CLb <<endl;
523  cout << "CLb inclusive [P(toy>=obs|0)] for observed 95% upper-limit "<< CLbinclusive <<endl;
524 
525  delete profile;
526  delete nll;
527 
528 }
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:777
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition: TSystem.cxx:1213
Bool_t saveSnapshot(const char *name, const char *paramNames)
Save snapshot of values and attributes (including "Constant") of parameters 'params' If importValue...
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:49
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:52
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:244
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found...
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:205
TCanvas * c1
Definition: legend1.C:2
ConfidenceBelt * GetConfidenceBelt()
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:265
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
std::vector< double > values
Definition: TwoHistoFit2D.C:32
const RooAbsCategoryLValue & indexCat() const
double confidenceLevel
static const char * filename()
#define gROOT
Definition: TROOT.h:344
virtual Double_t getMin(const char *name=0) const
int Int_t
Definition: RtypesCore.h:41
void OneSidedFrequentistUpperLimitWithBands(const char *infile="", const char *workspaceName="combined", const char *modelConfigName="ModelConfig", const char *dataName="obsData")
const TKDTreeBinning * bins
RooAbsPdf * getPdf(const char *catName) const
Return the p.d.f associated with the given index category name.
RooCmdArg Extended(Bool_t flag=kTRUE)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:250
Iterator abstract base class.
Definition: TIterator.h:32
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=1, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3851
RooAbsData * GetPointsToScan()
RooAbsArg * first() const
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:839
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:203
virtual void SetConfidenceLevel(Double_t cl)
set the confidence level for the interval (eg. 0.95 for a 95% Confidence Interval) ...
TText * tt
Definition: textangle.C:16
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2082
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
virtual void Draw(Option_t *option="")
Forward draw command to data store.
ConfidenceBelt is a concrete implementation of the ConfInterval interface.
RooCatType is an auxilary class for RooAbsCategory and defines a a single category state...
Definition: RooCatType.h:23
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
Double_t SignificanceToPValue(Double_t Z)
Definition: RooStatsUtils.h:64
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
tuple infile
Definition: mrt.py:15
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:203
virtual TestStatistic * GetTestStatistic(unsigned int i) const
Definition: ToyMCSampler.h:162
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
R__EXTERN TSystem * gSystem
Definition: TSystem.h:545
void SetNBins(Int_t bins)
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:293
void FluctuateNumDataEntries(bool flag=true)
void SetProofConfig(ProofConfig *pc=NULL)
Definition: ToyMCSampler.h:263
tuple w
Definition: qtexample.py:51
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
Bool_t loadSnapshot(const char *name)
Load the values and attributes of the parameters in the snapshot saved with the given name...
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:99
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:259
TestStatSampler * GetTestStatSampler() const
Returns instance of TestStatSampler.
virtual TObject * clone(const char *newname) const
Definition: RooArgSet.h:82
PointSetInterval is a concrete implementation of the ConfInterval interface.
tuple file
Definition: fildir.py:20
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:216
virtual PointSetInterval * GetInterval() const
Main interface to get a ConfInterval (will be a PointSetInterval)
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
virtual Double_t getMax(const char *name=0) const
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
virtual Double_t EvaluateTestStatistic(RooAbsData &data, RooArgSet &paramsOfInterest)=0
void Print(Option_t *opts=0) const
Print contents of the workspace.
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
Generate a new dataset containing the specified variables with events sampled from our distribution...
Definition: RooAbsPdf.cxx:1702
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:463
virtual TObject * Next()=0
#define NULL
Definition: Rtypes.h:82
void CreateConfBelt(bool flag=true)
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:247
Double_t GetAcceptanceRegionMax(RooArgSet &, Double_t cl=-1., Double_t leftside=-1.)
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event 'index'.
virtual const Text_t * GetName() const
Returns name of object.
Definition: RooCatType.h:45
TIterator * typeIterator() const
Return iterator over all defined states.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add element to non-owning set.
Definition: RooArgSet.cxx:448
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:42
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:527