DeanCLsExample/src/StandardHypoTestInvDemo.cc

// Standard tutorial macro for performing an inverted  hypothesis test 
//
// This macro will perform a scan of tehe p-values for computing the limit
// 

#include <iostream>

#include "StandardHypoTestInvDemo.h"


using namespace RooFit;
using namespace RooStats;


bool plotHypoTestResult = true; 
bool useProof = true;
bool optimize = false;
bool writeResult = false;
int nworkers = 4;
bool rebuild = false;
int nToyToRebuild = 100;


void StandardHypoTestInvDemo(const char * fileName,
                             const char * wsName,
                             const char * modelSBName,
                             const char * modelBName,
                             const char * dataName,
                             int calculatorType,
                             int testStatType, 
                             bool useCls,  
                             int npoints,   
                             double poimin,
                             double poimax,
                             int ntoys,
                             bool useNumberCounting,
                             const char * nuisPriorName)
{
/*

   Other Parameter to pass in tutorial
   apart from standard for filename, ws, modelconfig and data

    type = 0 Freq calculator 
    type = 1 Hybrid 

    testStatType = 0 LEP
                 = 1 Tevatron 
                 = 2 Profile Likelihood
                 = 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)

    useCLs          scan for CLs (otherwise for CLs+b)    

    npoints:        number of points to scan , for autoscan set npoints = -1 

    poimin,poimax:  min/max value to scan in case of fixed scans 
                    (if min >= max, try to find automatically)                           

    ntoys:         number of toys to use 

    useNumberCounting:  set to true when using number counting events 

    nuisPriorName:   name of prior for the nnuisance. This is often expressed as constraint term in the global model
                     It is needed only when using the HybridCalculator (type=1)
                     If not given by default the prior pdf from ModelConfig is used. 

    extra options are available as global paramters of the macro. They are: 

    plotHypoTestResult   plot result of tests at each point (TS distributions) 
    useProof = true;
    writeResult = true;
    nworkers = 4;


   */

   if (fileName==0) { 
      fileName = "example_combined_GaussExample_model.root";
      std::cout << "Use standard file generated with HistFactory :" << fileName << std::endl;
   }
   TFile * file = new TFile(fileName); 

   RooWorkspace * w = dynamic_cast<RooWorkspace*>( file->Get(wsName) );
   HypoTestInverterResult * r = 0; 
   std::cout << w << "\t" << fileName << std::endl;
   if (w != NULL) {
      r = RunInverter(w, modelSBName, modelBName, dataName, calculatorType, testStatType, npoints, poimin, poimax,  
                      ntoys, useCls, useNumberCounting, nuisPriorName );    
      if (!r) { 
         std::cerr << "Error running the HypoTestInverter - Exit " << std::endl;
         return;          
      }
   }
   else 
   { 
      // case workspace is not present look for the inverter result
      std::cout << "Reading an HypoTestInverterResult with name " << wsName << " from file " << fileName << std::endl;
      r = dynamic_cast<HypoTestInverterResult*>( file->Get(wsName) ); //
      if (!r) { 
         std::cerr << "File " << fileName << " does not contain a workspace or an HypoTestInverterResult - Exit " 
                   << std::endl;
         file->ls();
         return; 
      }
   }    
          

   double upperLimit = r->UpperLimit();
   double ulError = r->UpperLimitEstimatedError();


   std::cout << "The computed upper limit is: " << upperLimit << " +/- " << ulError << std::endl;
 
   const int nEntries = r->ArraySize();


   const char *  typeName = (calculatorType == 0) ? "Frequentist" : "Hybrid";
   const char * resultName = (w) ? w->GetName() : r->GetName();
   TString plotTitle = TString::Format("%s CL Scan for workspace %s",typeName,resultName);
   HypoTestInverterPlot *plot = new HypoTestInverterPlot("HTI_Result_Plot",plotTitle,r);
   plot->Draw("CLb 2CL");  // plot all and Clb

   if (plotHypoTestResult) { 
      TCanvas * c2 = new TCanvas();
      c2->Divide( 2, TMath::Ceil(nEntries/2));
      for (int i=0; i<nEntries; i++) {
         c2->cd(i+1);
         SamplingDistPlot * pl = plot->MakeTestStatPlot(i);
         pl->SetLogYaxis(true);
         pl->Draw();
      }
   }


   std::cout << " expected limit (median) " <<  r->GetExpectedUpperLimit(0) << std::endl;
   std::cout << " expected limit (-1 sig) " << r->GetExpectedUpperLimit(-1) << std::endl;
   std::cout << " expected limit (+1 sig) " << r->GetExpectedUpperLimit(1) << std::endl;


   if (w != NULL && writeResult) {

      // write to a file the results
      const char *  calcType = (calculatorType == 0) ? "Freq" : "Hybr";
      const char *  limitType = (useCls) ? "CLs" : "Cls+b";
      const char * scanType = (npoints < 0) ? "auto" : "grid";
      TString resultFileName = TString::Format("%s_%s_%s_ts%d_",calcType,limitType,scanType,testStatType);      
      resultFileName += fileName;
      
      TFile * fileOut = new TFile(resultFileName,"RECREATE");
      r->Write();
      fileOut->Close();                                                                     
   }   

}


// internal routine to run the inverter
HypoTestInverterResult *  RunInverter(RooWorkspace * w, const char * modelSBName, const char * modelBName, 
                                      const char * dataName, int type,  int testStatType, 
                                      int npoints, double poimin, double poimax, 
                                      int ntoys, bool useCls, bool useNumberCounting, const char * nuisPriorName ) 
{

   std::cout << "Running HypoTestInverter on the workspace " << w->GetName() << std::endl;

   w->Print();


   RooAbsData * data = w->data(dataName); 
   if (!data) { 
      Error("StandardHypoTestDemo","Not existing data %s",dataName);
      return 0;
   }
   else 
      std::cout << "Using data set " << dataName << std::endl;

   
   // get models from WS
   // get the modelConfig out of the file
   ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);
   ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);

   if (!sbModel) {
      Error("StandardHypoTestDemo","Not existing ModelConfig %s",modelSBName);
      return 0;
   }
   // check the model 
   if (!sbModel->GetPdf()) { 
      Error("StandardHypoTestDemo","Model %s has no pdf ",modelSBName);
      return 0;
   }
   if (!sbModel->GetParametersOfInterest()) {
      Error("StandardHypoTestDemo","Model %s has no poi ",modelSBName);
      return 0;
   }
   if (!sbModel->GetParametersOfInterest()) {
      Error("StandardHypoTestInvDemo","Model %s has no poi ",modelSBName);
      return 0;
   }
   if (!sbModel->GetSnapshot() ) { 
      Info("StandardHypoTestInvDemo","Model %s has no snapshot  - make one using model poi",modelSBName);
      sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );
   }


   if (!bModel || bModel == sbModel) {
      Info("StandardHypoTestInvDemo","The background model %s does not exist",modelBName);
      Info("StandardHypoTestInvDemo","Copy it from ModelConfig %s and set POI to zero",modelSBName);
      bModel = (ModelConfig*) sbModel->Clone();
      bModel->SetName(TString(modelSBName)+TString("_with_poi_0"));      
      RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
      if (!var) return 0;
      double oldval = var->getVal();
      var->setVal(0);
      bModel->SetSnapshot( RooArgSet(*var)  );
      var->setVal(oldval);
   }
   else { 
      if (!bModel->GetSnapshot() ) { 
         Info("StandardHypoTestInvDemo","Model %s has no snapshot  - make one using model poi and 0 values ",modelBName);
         RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
         if (var) { 
            double oldval = var->getVal();
            var->setVal(0);
            bModel->SetSnapshot( RooArgSet(*var)  );
            var->setVal(oldval);
         }
         else { 
            Error("StandardHypoTestInvDemo","Model %s has no valid poi",modelBName);
            return 0;
         }         
      }
   }


   SimpleLikelihoodRatioTestStat slrts(*sbModel->GetPdf(),*bModel->GetPdf());
   if (sbModel->GetSnapshot()) slrts.SetNullParameters(*sbModel->GetSnapshot());
   if (bModel->GetSnapshot()) slrts.SetAltParameters(*bModel->GetSnapshot());

   // ratio of profile likelihood - need to pass snapshot for the alt
   RatioOfProfiledLikelihoodsTestStat 
      ropl(*sbModel->GetPdf(), *bModel->GetPdf(), bModel->GetSnapshot());
   ropl.SetSubtractMLE(false);
   
   ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   if (testStatType == 3) profll.SetOneSided(1);
   if (optimize) profll.SetReuseNLL(true);

   TestStatistic * testStat = &slrts;
   if (testStatType == 1) testStat = &ropl;
   if (testStatType == 2 || testStatType == 3) testStat = &profll;
  
   
   HypoTestCalculatorGeneric *  hc = 0;
   if (type == 0) hc = new FrequentistCalculator(*data, *bModel, *sbModel);
   else hc = new HybridCalculator(*data, *bModel, *sbModel);

   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   if (useNumberCounting) toymcs->SetNEventsPerToy(1);
   toymcs->SetTestStatistic(testStat);
   if (optimize) toymcs->SetUseMultiGen(true);


   if (type == 1) { 
      HybridCalculator *hhc = (HybridCalculator*) hc;
      hhc->SetToys(ntoys,ntoys); 

      // check for nuisance prior pdf in case of nuisance parameters 
      if (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters() ) {
         RooAbsPdf * nuisPdf = 0; 
         if (nuisPriorName) nuisPdf = w->pdf(nuisPriorName);
         // use prior defined first in bModel (then in SbModel)
         if (!nuisPdf)  { 
            Info("StandardHypoTestInvDemo","No nuisance pdf given for the HybridCalculator - try to use the prior pdf from the model");
            nuisPdf = (bModel->GetPriorPdf() ) ?  bModel->GetPriorPdf() : sbModel->GetPriorPdf();
         }
         if (!nuisPdf) { 
            Error("StandardHypoTestInvDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified");
            return 0;
         }
         const RooArgSet * nuisParams = (bModel->GetNuisanceParameters() ) ? bModel->GetNuisanceParameters() : sbModel->GetNuisanceParameters();
         RooArgSet * np = nuisPdf->getObservables(*nuisParams);
         if (np->getSize() == 0) { 
            Warning("StandardHypoTestInvDemo","Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range");
         }
         delete np;
         hhc->ForcePriorNuisanceAlt(*nuisPdf);
         hhc->ForcePriorNuisanceNull(*nuisPdf);
      }
   } 
   else 
      ((FrequentistCalculator*) hc)->SetToys(ntoys,ntoys); 

   // Get the result
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);


   TStopwatch tw; tw.Start(); 
   const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
   RooRealVar *poi = (RooRealVar*)poiSet->first();

   // fit the data first
   sbModel->GetPdf()->fitTo(*data);
   double poihat  = poi->getVal();


   HypoTestInverter calc(*hc);
   calc.SetConfidenceLevel(0.95);

   calc.UseCLs(useCls);
   calc.SetVerbose(true);

   // can speed up using proof-lite
   if (useProof && nworkers > 1) { 
      ProofConfig pc(*w, nworkers, "", kFALSE);
      toymcs->SetProofConfig(&pc);    // enable proof
   }

   
   if (npoints > 0) {
      if (poimin >= poimax) { 
         // if no min/max given scan between MLE and +4 sigma 
         poimin = int(poihat);
         poimax = int(poihat +  4 * poi->getError());
      }
      std::cout << "Doing a fixed scan  in interval : " << poimin << " , " << poimax << std::endl;
      //calc.SetFixedScan(1,6,6);
      calc.SetFixedScan(npoints,poimin,poimax);
   }
   else { 
      //poi->setMax(10*int( (poihat+ 10 *poi->getError() )/10 ) );
      std::cout << "Doing an  automatic scan  in interval : " << poi->getMin() << " , " << poi->getMax() << std::endl;
   }

   HypoTestInverterResult * r = calc.GetInterval();

   if (rebuild) {
      calc.SetCloseProof(1);
      SamplingDistribution * limDist = calc.GetUpperLimitDistribution(true,nToyToRebuild);
      if (limDist) { 
         std::cout << "expected up limit " << limDist->InverseCDF(0.5) << " +/- " 
                   << limDist->InverseCDF(0.16) << "  " 
                   << limDist->InverseCDF(0.84) << "\n"; 
      }
      else 
         std::cout << "ERROR : failed to re-build distributions " << std::endl; 
   }

   //update r to a new re-freshed copied
   r = calc.GetInterval();

   delete hc;

   return r; 
}

void ReadResult(const char * fileName, const char * resultName="", bool useCLs=true) { 
   // read a previous stored result from a file given the result name

   StandardHypoTestInvDemo(fileName, resultName,"","","",0,0,useCLs);
}

Revision:	1.1
Committed:	Thu Sep 29 14:52:09 2011 UTC (13 years, 7 months ago) by dhidas
Content type:	text/plain
Branch:	MAIN
Branch point for:	dhidas
Log Message:	Initial revision
#	Content
1	// Standard tutorial macro for performing an inverted hypothesis test
2	//
3	// This macro will perform a scan of tehe p-values for computing the limit
4	//
5
6	#include <iostream>
7
8	#include "StandardHypoTestInvDemo.h"
9
10
11	using namespace RooFit;
12	using namespace RooStats;
13
14
15	bool plotHypoTestResult = true;
16	bool useProof = true;
17	bool optimize = false;
18	bool writeResult = false;
19	int nworkers = 4;
20	bool rebuild = false;
21	int nToyToRebuild = 100;
22
23
24
25
26
27	void StandardHypoTestInvDemo(const char * fileName,
28	const char * wsName,
29	const char * modelSBName,
30	const char * modelBName,
31	const char * dataName,
32	int calculatorType,
33	int testStatType,
34	bool useCls,
35	int npoints,
36	double poimin,
37	double poimax,
38	int ntoys,
39	bool useNumberCounting,
40	const char * nuisPriorName)
41	{
42	/*
43
44	Other Parameter to pass in tutorial
45	apart from standard for filename, ws, modelconfig and data
46
47	type = 0 Freq calculator
48	type = 1 Hybrid
49
50	testStatType = 0 LEP
51	= 1 Tevatron
52	= 2 Profile Likelihood
53	= 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)
54
55	useCLs scan for CLs (otherwise for CLs+b)
56
57	npoints: number of points to scan , for autoscan set npoints = -1
58
59	poimin,poimax: min/max value to scan in case of fixed scans
60	(if min >= max, try to find automatically)
61
62	ntoys: number of toys to use
63
64	useNumberCounting: set to true when using number counting events
65
66	nuisPriorName: name of prior for the nnuisance. This is often expressed as constraint term in the global model
67	It is needed only when using the HybridCalculator (type=1)
68	If not given by default the prior pdf from ModelConfig is used.
69
70	extra options are available as global paramters of the macro. They are:
71
72	plotHypoTestResult plot result of tests at each point (TS distributions)
73	useProof = true;
74	writeResult = true;
75	nworkers = 4;
76
77
78	*/
79
80	if (fileName==0) {
81	fileName = "example_combined_GaussExample_model.root";
82	std::cout << "Use standard file generated with HistFactory :" << fileName << std::endl;
83	}
84	TFile * file = new TFile(fileName);
85
86	RooWorkspace * w = dynamic_cast<RooWorkspace*>( file->Get(wsName) );
87	HypoTestInverterResult * r = 0;
88	std::cout << w << "\t" << fileName << std::endl;
89	if (w != NULL) {
90	r = RunInverter(w, modelSBName, modelBName, dataName, calculatorType, testStatType, npoints, poimin, poimax,
91	ntoys, useCls, useNumberCounting, nuisPriorName );
92	if (!r) {
93	std::cerr << "Error running the HypoTestInverter - Exit " << std::endl;
94	return;
95	}
96	}
97	else
98	{
99	// case workspace is not present look for the inverter result
100	std::cout << "Reading an HypoTestInverterResult with name " << wsName << " from file " << fileName << std::endl;
101	r = dynamic_cast<HypoTestInverterResult*>( file->Get(wsName) ); //
102	if (!r) {
103	std::cerr << "File " << fileName << " does not contain a workspace or an HypoTestInverterResult - Exit "
104	<< std::endl;
105	file->ls();
106	return;
107	}
108	}
109
110
111	double upperLimit = r->UpperLimit();
112	double ulError = r->UpperLimitEstimatedError();
113
114
115	std::cout << "The computed upper limit is: " << upperLimit << " +/- " << ulError << std::endl;
116
117	const int nEntries = r->ArraySize();
118
119
120	const char * typeName = (calculatorType == 0) ? "Frequentist" : "Hybrid";
121	const char * resultName = (w) ? w->GetName() : r->GetName();
122	TString plotTitle = TString::Format("%s CL Scan for workspace %s",typeName,resultName);
123	HypoTestInverterPlot *plot = new HypoTestInverterPlot("HTI_Result_Plot",plotTitle,r);
124	plot->Draw("CLb 2CL"); // plot all and Clb
125
126	if (plotHypoTestResult) {
127	TCanvas * c2 = new TCanvas();
128	c2->Divide( 2, TMath::Ceil(nEntries/2));
129	for (int i=0; i<nEntries; i++) {
130	c2->cd(i+1);
131	SamplingDistPlot * pl = plot->MakeTestStatPlot(i);
132	pl->SetLogYaxis(true);
133	pl->Draw();
134	}
135	}
136
137
138	std::cout << " expected limit (median) " << r->GetExpectedUpperLimit(0) << std::endl;
139	std::cout << " expected limit (-1 sig) " << r->GetExpectedUpperLimit(-1) << std::endl;
140	std::cout << " expected limit (+1 sig) " << r->GetExpectedUpperLimit(1) << std::endl;
141
142
143	if (w != NULL && writeResult) {
144
145	// write to a file the results
146	const char * calcType = (calculatorType == 0) ? "Freq" : "Hybr";
147	const char * limitType = (useCls) ? "CLs" : "Cls+b";
148	const char * scanType = (npoints < 0) ? "auto" : "grid";
149	TString resultFileName = TString::Format("%s_%s_%s_ts%d_",calcType,limitType,scanType,testStatType);
150	resultFileName += fileName;
151
152	TFile * fileOut = new TFile(resultFileName,"RECREATE");
153	r->Write();
154	fileOut->Close();
155	}
156
157	}
158
159
160	// internal routine to run the inverter
161	HypoTestInverterResult * RunInverter(RooWorkspace * w, const char * modelSBName, const char * modelBName,
162	const char * dataName, int type, int testStatType,
163	int npoints, double poimin, double poimax,
164	int ntoys, bool useCls, bool useNumberCounting, const char * nuisPriorName )
165	{
166
167	std::cout << "Running HypoTestInverter on the workspace " << w->GetName() << std::endl;
168
169	w->Print();
170
171
172	RooAbsData * data = w->data(dataName);
173	if (!data) {
174	Error("StandardHypoTestDemo","Not existing data %s",dataName);
175	return 0;
176	}
177	else
178	std::cout << "Using data set " << dataName << std::endl;
179
180
181	// get models from WS
182	// get the modelConfig out of the file
183	ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);
184	ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);
185
186	if (!sbModel) {
187	Error("StandardHypoTestDemo","Not existing ModelConfig %s",modelSBName);
188	return 0;
189	}
190	// check the model
191	if (!sbModel->GetPdf()) {
192	Error("StandardHypoTestDemo","Model %s has no pdf ",modelSBName);
193	return 0;
194	}
195	if (!sbModel->GetParametersOfInterest()) {
196	Error("StandardHypoTestDemo","Model %s has no poi ",modelSBName);
197	return 0;
198	}
199	if (!sbModel->GetParametersOfInterest()) {
200	Error("StandardHypoTestInvDemo","Model %s has no poi ",modelSBName);
201	return 0;
202	}
203	if (!sbModel->GetSnapshot() ) {
204	Info("StandardHypoTestInvDemo","Model %s has no snapshot - make one using model poi",modelSBName);
205	sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );
206	}
207
208
209	if (!bModel \|\| bModel == sbModel) {
210	Info("StandardHypoTestInvDemo","The background model %s does not exist",modelBName);
211	Info("StandardHypoTestInvDemo","Copy it from ModelConfig %s and set POI to zero",modelSBName);
212	bModel = (ModelConfig*) sbModel->Clone();
213	bModel->SetName(TString(modelSBName)+TString("_with_poi_0"));
214	RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
215	if (!var) return 0;
216	double oldval = var->getVal();
217	var->setVal(0);
218	bModel->SetSnapshot( RooArgSet(*var) );
219	var->setVal(oldval);
220	}
221	else {
222	if (!bModel->GetSnapshot() ) {
223	Info("StandardHypoTestInvDemo","Model %s has no snapshot - make one using model poi and 0 values ",modelBName);
224	RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
225	if (var) {
226	double oldval = var->getVal();
227	var->setVal(0);
228	bModel->SetSnapshot( RooArgSet(*var) );
229	var->setVal(oldval);
230	}
231	else {
232	Error("StandardHypoTestInvDemo","Model %s has no valid poi",modelBName);
233	return 0;
234	}
235	}
236	}
237
238
239	SimpleLikelihoodRatioTestStat slrts(sbModel->GetPdf(),bModel->GetPdf());
240	if (sbModel->GetSnapshot()) slrts.SetNullParameters(*sbModel->GetSnapshot());
241	if (bModel->GetSnapshot()) slrts.SetAltParameters(*bModel->GetSnapshot());
242
243	// ratio of profile likelihood - need to pass snapshot for the alt
244	RatioOfProfiledLikelihoodsTestStat
245	ropl(sbModel->GetPdf(), bModel->GetPdf(), bModel->GetSnapshot());
246	ropl.SetSubtractMLE(false);
247
248	ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
249	if (testStatType == 3) profll.SetOneSided(1);
250	if (optimize) profll.SetReuseNLL(true);
251
252	TestStatistic * testStat = &slrts;
253	if (testStatType == 1) testStat = &ropl;
254	if (testStatType == 2 \|\| testStatType == 3) testStat = &profll;
255
256
257	HypoTestCalculatorGeneric * hc = 0;
258	if (type == 0) hc = new FrequentistCalculator(data, bModel, *sbModel);
259	else hc = new HybridCalculator(data, bModel, *sbModel);
260
261	ToyMCSampler toymcs = (ToyMCSampler)hc->GetTestStatSampler();
262	if (useNumberCounting) toymcs->SetNEventsPerToy(1);
263	toymcs->SetTestStatistic(testStat);
264	if (optimize) toymcs->SetUseMultiGen(true);
265
266
267	if (type == 1) {
268	HybridCalculator hhc = (HybridCalculator) hc;
269	hhc->SetToys(ntoys,ntoys);
270
271	// check for nuisance prior pdf in case of nuisance parameters
272	if (bModel->GetNuisanceParameters() \|\| sbModel->GetNuisanceParameters() ) {
273	RooAbsPdf * nuisPdf = 0;
274	if (nuisPriorName) nuisPdf = w->pdf(nuisPriorName);
275	// use prior defined first in bModel (then in SbModel)
276	if (!nuisPdf) {
277	Info("StandardHypoTestInvDemo","No nuisance pdf given for the HybridCalculator - try to use the prior pdf from the model");
278	nuisPdf = (bModel->GetPriorPdf() ) ? bModel->GetPriorPdf() : sbModel->GetPriorPdf();
279	}
280	if (!nuisPdf) {
281	Error("StandardHypoTestInvDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified");
282	return 0;
283	}
284	const RooArgSet * nuisParams = (bModel->GetNuisanceParameters() ) ? bModel->GetNuisanceParameters() : sbModel->GetNuisanceParameters();
285	RooArgSet * np = nuisPdf->getObservables(*nuisParams);
286	if (np->getSize() == 0) {
287	Warning("StandardHypoTestInvDemo","Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range");
288	}
289	delete np;
290	hhc->ForcePriorNuisanceAlt(*nuisPdf);
291	hhc->ForcePriorNuisanceNull(*nuisPdf);
292	}
293	}
294	else
295	((FrequentistCalculator*) hc)->SetToys(ntoys,ntoys);
296
297	// Get the result
298	RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
299
300
301	TStopwatch tw; tw.Start();
302	const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
303	RooRealVar poi = (RooRealVar)poiSet->first();
304
305	// fit the data first
306	sbModel->GetPdf()->fitTo(*data);
307	double poihat = poi->getVal();
308
309
310	HypoTestInverter calc(*hc);
311	calc.SetConfidenceLevel(0.95);
312
313	calc.UseCLs(useCls);
314	calc.SetVerbose(true);
315
316	// can speed up using proof-lite
317	if (useProof && nworkers > 1) {
318	ProofConfig pc(*w, nworkers, "", kFALSE);
319	toymcs->SetProofConfig(&pc); // enable proof
320	}
321
322
323	if (npoints > 0) {
324	if (poimin >= poimax) {
325	// if no min/max given scan between MLE and +4 sigma
326	poimin = int(poihat);
327	poimax = int(poihat + 4 * poi->getError());
328	}
329	std::cout << "Doing a fixed scan in interval : " << poimin << " , " << poimax << std::endl;
330	//calc.SetFixedScan(1,6,6);
331	calc.SetFixedScan(npoints,poimin,poimax);
332	}
333	else {
334	//poi->setMax(10int( (poihat+ 10 poi->getError() )/10 ) );
335	std::cout << "Doing an automatic scan in interval : " << poi->getMin() << " , " << poi->getMax() << std::endl;
336	}
337
338	HypoTestInverterResult * r = calc.GetInterval();
339
340	if (rebuild) {
341	calc.SetCloseProof(1);
342	SamplingDistribution * limDist = calc.GetUpperLimitDistribution(true,nToyToRebuild);
343	if (limDist) {
344	std::cout << "expected up limit " << limDist->InverseCDF(0.5) << " +/- "
345	<< limDist->InverseCDF(0.16) << " "
346	<< limDist->InverseCDF(0.84) << "\n";
347	}
348	else
349	std::cout << "ERROR : failed to re-build distributions " << std::endl;
350	}
351
352	//update r to a new re-freshed copied
353	r = calc.GetInterval();
354
355	delete hc;
356
357	return r;
358	}
359
360	void ReadResult(const char * fileName, const char * resultName="", bool useCLs=true) {
361	// read a previous stored result from a file given the result name
362
363	StandardHypoTestInvDemo(fileName, resultName,"","","",0,0,useCLs);
364	}
365