C++ RooFitResult类代码示例

void voigtian(RooDataSet *dataset, RooDataSet *dataset2, RooRealVar &variable, RooPlot *fitFrame, RooBinning b, double rangeMin, double rangeMax, vector <double> &fitParameters)
{
        RooRealVar mean("mean","mean",0.0,-0.1,0.1);
        RooRealVar sigma("sigma","sigma",0.5,0.0,1.0);
        RooRealVar width("width","width",0.5,0.0,1.0);

        RooVoigtian * pdf = new RooVoigtian("pdf","Voigtian",variable,mean,sigma,width);

        dataset->plotOn(fitFrame,Name("myhist"),Binning(b),DataError(RooAbsData::SumW2));

        RooFitResult * res = pdf->fitTo(*dataset, Range(rangeMin, rangeMax),Save(),SumW2Error(kTRUE));
	res->Print();

	//minNll = res->minNll();
        
	pdf->plotOn(fitFrame,Name("mycurve"));


	fitParameters.push_back(3); // nb of fit parameters

	fitParameters.push_back(mean.getVal());
	fitParameters.push_back(mean.getError());

	fitParameters.push_back(sigma.getVal());
        fitParameters.push_back(sigma.getError());
	
	fitParameters.push_back(width.getVal());
        fitParameters.push_back(width.getError());		
	

}

void PDF_GLWADS_Dpi_K3pi::setUncertainties(config c)
{
    switch(c)
    {
    case lumi1fb:
    {
        obsErrSource = "1fb-1, ExpNll/sept2012K3PIResult.root";
        TString File = this->dir+"/ExpNll/sept2012K3PIResult.root";
        TFile *fr = TFile::Open(File);
        RooFitResult *r = (RooFitResult*)fr->Get("fitresult_model_reducedData_binned");
        assert(r);
        for ( int i=0; i<nObs; i++ )
        {
            RooRealVar* pObs = (RooRealVar*)((RooArgList*)observables)->at(i);
            RooRealVar* pRes = (RooRealVar*)r->floatParsFinal().find(obsTmkToMalcolm(pObs->GetName()));
            assert(pRes);
            StatErr[i] = pRes->getError();
        }
        SystErr[0] = 0.010;   // afav_dpi_obs
        SystErr[1] = 0.00011; // rp_dpi_obs
        SystErr[2] = 0.00011; // rm_dpi_obs
        fr->Close();
        delete r;
        delete fr;
        break;
    }
    default:
        cout << "PDF_GLWADS_Dpi_K3pi::setUncertainties() : ERROR : config "+ConfigToTString(c)+" not found." << endl;
        exit(1);
    }
}

RooFitResult *breakDownFit(RooSimultaneous *m, RooAbsData *d, RooRealVar *mass, bool precondition = false){
	if(precondition){
		 const char *catsName = m->indexCat().GetName();
		 TIterator *it = m->indexCat().typeIterator();
		 while(RooCatType* ci = dynamic_cast<RooCatType*>(it->Next())) {
			 const Text_t *catLabel = ci->GetName();
			 RooAbsPdf *pdf = m->getPdf(Form("%s",catLabel));
			 RooAbsData *reduced = d->reduce(SelectVars(*mass),Cut(Form("%s==%s::%s",catsName, catsName, catLabel)));
			 RooFitResult *r = pdf->fitTo(*reduced,PrintLevel(-1),Save(),
					 Minimizer("Minuit2","migrad"),Strategy(0),Hesse(false),Minos(false),Optimize(false)
					 );
			 cout << catsName << " " << catLabel << " M2migrad0 " << r->status() << endl;
			 if(r->status()!=0){
				 RooFitResult *r = pdf->fitTo(*reduced, PrintLevel(-1), Save());
				 cout << catsName << " " << catLabel << " Mmigrad1 " << r->status() << endl;
			 }
		 }
	}

	RooFitResult *r = m->fitTo(*d, Save(), PrintLevel(-1),
			Strategy(0));
	cout << "Global fit Mmigrad0 " << r->status() << endl;
	if(r->status()!=0){
	 RooFitResult *r = m->fitTo(*d, PrintLevel(-1), Save(),
			 Minimizer("Minuit","minimize"),Strategy(2));
	 cout << "Global fit Mminimize2 " << r->status() << endl;
	 return r;
	}

	return r;
}

void PDF_GLWADS_DKDpi_K3pi::setObservables(config c)
{
	switch(c)
	{
		case truth:{
					   setObservablesTruth();
					   break;
				   }
		case toy:{
					 setObservablesToy();
					 break;
				 }
		case lumi1fb:{
						 obsValSource = "1fb-1, ExpNll/sept2012K3PIResult.root";
						 TString File = this->dir+"/ExpNll/sept2012K3PIResult.root";
						 TFile *fr = TFile::Open(File);
						 RooFitResult *r = (RooFitResult*)fr->Get("fitresult_model_reducedData_binned");
						 assert(r);
						 TIterator* it = observables->createIterator();
						 while ( RooRealVar* pObs = (RooRealVar*)it->Next() )
						 {
							 RooRealVar* pRes = (RooRealVar*)r->floatParsFinal().find(obsTmkToMalcolm(pObs->GetName()));
							 pObs->setVal(pRes->getVal());
						 }
						 fr->Close();
						 delete r;
						 delete fr;
						 break;
					 }
		case lumi3fb:{
						 obsValSource = "3fb-1 ANA v7 unblind"; // https://twiki.cern.ch/twiki/pub/LHCbPhysics/B2D0K/LHCb-ANA-2014-071-v7.pdf (see Vavas email 04/08/15)

						 // these get transformed over from the new inputs using ExpNll/transportGLWADS_new_to_old.py
						 // in the case of the DK only (robust) combination some of the observables don't exist
						 // usemap as the temp store
						 std::map< TString, double > vals;
						 vals["rkp_k3pi_obs"]      =  0.0793;
						 vals["afav_dk_k3pi_obs"]  =  -0.0004;
						 vals["afav_dpi_k3pi_obs"] =  0.0;
						 vals["rp_dk_k3pi_obs"]    =  0.018369;
						 vals["rm_dk_k3pi_obs"]    =  0.009611;
						 vals["rp_dpi_k3pi_obs"]   =  0.003683;
						 vals["rm_dpi_k3pi_obs"]   =  0.003857;

						 // now can loop the observables and set the values
						 TIterator* it = observables->createIterator();
						 while ( RooRealVar* pObs = (RooRealVar*)it->Next() ){
							 pObs->setVal(vals[pObs->GetName()]);
						 }

						 vals.clear();
						 break;
				}
		default:{
					cout << "PDF_GLWADS_DKDpi_K3pi::setObservables() : ERROR : config "+ConfigToTString(c)+" not found." << endl;
					exit(1);
				}
	}
}

void checkBestFitPoint(std::string workspace, std::string fitFile, bool splusb){
	
	// Open the ws file...
	TFile *fd_=0;
	TFile *fw_=0;
		
	gSystem->Load("$CMSSW_BASE/lib/$SCRAM_ARCH/libHiggsAnalysisCombinedLimit.so");
	gROOT->SetBatch(true);
	gStyle->SetOptFit(0);
	gStyle->SetOptStat(0);
	gStyle->SetPalette(1,0);

	fw_ =  TFile::Open(workspace.c_str());
	w   = (RooWorkspace*) fw_->Get("w");
  w->Print();
	RooDataSet *data = (RooDataSet*) w->data("data_obs");
  if (splusb) {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig");
  } else {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig_bonly");
  }
	std::cout << "make nll"<<std::endl;
	nll = mc_s->GetPdf()->createNLL(
		*data,RooFit::Constrain(*mc_s->GetNuisanceParameters())
		,RooFit::Extended(mc_s->GetPdf()->canBeExtended()));
	
	// Now get the best fit result
	fd_ =  TFile::Open(fitFile.c_str());
	RooFitResult *fit;
  if (splusb) {
    fit =(RooFitResult*)fd_->Get("fit_s");
  } else {
    fit =(RooFitResult*)fd_->Get("fit_b");
  }
	RooArgSet fitargs = fit->floatParsFinal();
	
	std::cout << "Got the best fit values" <<std::endl;		
	w->saveSnapshot("bestfitall",fitargs,true);
	
  TString filename;
  if (splusb) {
    filename = "minimum_s.pdf";
  } else {
    filename = "minimum_b.pdf";
  }
	// Now make the plots!	
	TCanvas *c = new TCanvas("c","",600,600);
	c->SaveAs((filename+"["));

	TIterator* iter(fitargs->createIterator());
        for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
                 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);      
                 std::string name = rrv->GetName();
		 TGraph *gr = graphLH(name,rrv->getError());
		 gr->Draw("ALP");
		 c->SaveAs((filename+"["));
	}
	c->SaveAs((filename+"]"));
}

float getDLL(RooWorkspace* w, TString tag) {
  //RooFitResult::sexp_b_fitres

  RooFitResult *b = (RooFitResult*)w->obj(tag+"_b_fitres");
  RooFitResult *sb = (RooFitResult*)w->obj(tag+"_sb_fitres");

  return -2*(sb->minNll()-b->minNll());
}

void rf314_paramfitrange()
{

  // D e f i n e   o b s e r v a b l e s   a n d   d e c a y   p d f 
  // ---------------------------------------------------------------

  // Declare observables
  RooRealVar t("t","t",0,5) ;
  RooRealVar tmin("tmin","tmin",0,0,5) ;

  // Make parameterized range in t : [tmin,5]
  t.setRange(tmin,RooConst(t.getMax())) ;

  // Make pdf
  RooRealVar tau("tau","tau",-1.54,-10,-0.1) ;
  RooExponential model("model","model",t,tau) ;



  // C r e a t e   i n p u t   d a t a 
  // ------------------------------------

  // Generate complete dataset without acceptance cuts (for reference)
  RooDataSet* dall = model.generate(t,10000) ;

  // Generate a (fake) prototype dataset for acceptance limit values
  RooDataSet* tmp = RooGaussian("gmin","gmin",tmin,RooConst(0),RooConst(0.5)).generate(tmin,5000) ;

  // Generate dataset with t values that observe (t>tmin)
  RooDataSet* dacc = model.generate(t,ProtoData(*tmp)) ;



  // F i t   p d f   t o   d a t a   i n   a c c e p t a n c e   r e g i o n
  // -----------------------------------------------------------------------

  RooFitResult* r = model.fitTo(*dacc,Save()) ;


 
  // P l o t   f i t t e d   p d f   o n   f u l l   a n d   a c c e p t e d   d a t a 
  // ---------------------------------------------------------------------------------

  // Make plot frame, add datasets and overlay model
  RooPlot* frame = t.frame(Title("Fit to data with per-event acceptance")) ;
  dall->plotOn(frame,MarkerColor(kRed),LineColor(kRed)) ;
  model.plotOn(frame) ;
  dacc->plotOn(frame) ;

  // Print fit results to demonstrate absence of bias
  r->Print("v") ;


  new TCanvas("rf314_paramranges","rf314_paramranges",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;

  return ;
}

MCMCInterval * Tprime::GetMcmcInterval(ModelConfig mc,
                                       double conf_level,
                                       int n_iter,
                                       int n_burn,
                                       double left_side_tail_fraction,
                                       int n_bins) {
    //
    // Bayesian MCMC calculation using arbitrary ModelConfig
    // Want an efficient proposal function, so derive it from covariance
    // matrix of fit
    //

    RooAbsData * _data = data;
    //RooAbsData * _data = pWs->data("obsData");
    //RooStats::ModelConfig * _mc = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");
    RooStats::ModelConfig * _mc = GetModelConfig();
    _mc->Print();

    //RooFitResult * fit = pWs->pdf("model_tprime")->fitTo(*_data,Save());
    RooFitResult * fit = _mc->GetPdf()->fitTo(*_data,Save());
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
    ph.SetCacheSize(100);
    ProposalFunction * pf = ph.GetProposalFunction();

    //delete pf;
    //pf = new SequentialProposal();

    MCMCCalculator mcmc( *_data, mc );
    mcmc.SetConfidenceLevel(conf_level);
    mcmc.SetNumIters(n_iter);          // Metropolis-Hastings algorithm iterations
    mcmc.SetProposalFunction(*pf);
    mcmc.SetNumBurnInSteps(n_burn); // first N steps to be ignored as burn-in
    mcmc.SetLeftSideTailFraction(left_side_tail_fraction);
    mcmc.SetNumBins(n_bins);

    //mcInt = mcmc.GetInterval();
    try {
        mcInt = mcmc.GetInterval();
    } catch ( std::length_error &ex) {
        mcInt = 0;
    }

    //std::cout << "!!!!!!!!!!!!!! interval" << std::endl;
    if (mcInt == 0) std::cout << "No interval found!" << std::endl;

    delete fit;
    delete pf;

    return mcInt;
}

//_________________________________________________
void TestJeffreysGaussSigma(){
  // this one is VERY sensitive
  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
  //   and you get really bizzare shapes
  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
  //   and the PDF falls off too fast at high sigma
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
  w.factory("n[100,.1,2000]");
  w.factory("ExtendPdf::p(g,n)");
  //  w.var("sigma")->setConstant();
  w.var("mu")->setConstant();
  w.var("n")->setConstant();
  w.var("x")->setBins(301);

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;


  //  w.defineSet("poi","mu,sigma");
  //w.defineSet("poi","mu,sigma,n");
  w.defineSet("poi","sigma");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();
  //  return;

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","sqrt(2.)/sigma",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("sigma")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

void fit(Int_t i, Double_t va=-0.43, Double_t vb=0.){
    sprintf(namestr,"%d_%.2f_%.2f",i,va,vb);

    Double_t g1, g2;
    g1 = 0.1*r.Rndm() - 0.05;
    g2 =     r.Rndm() - 0.5;

    //a->setVal(va);
    //b->setVal(vb);

    //G001->setVal(g1);
    //G002->setVal(g2);

    //RooRealVar * G000 = new RooRealVar("G000", "G000",  0.5);
    //RooRealVar * G001 = new RooRealVar("G001", "G001",  g1, -5., 5.);
    //RooRealVar * G002 = new RooRealVar("G002", "G002",  g2, -5., 5.);
    //RooRealVar * G003 = new RooRealVar("G003", "G003", 0.0);//, -1., 1.);
    //RooRealVar * G004 = new RooRealVar("G004", "G004", 0.0);//, -1., 1.);
    //RooRealVar * a    = new RooRealVar("a",    "a",    va);
    //RooRealVar * b    = new RooRealVar("b",    "b",    vb);
    //RooRealVar * n    = new RooRealVar("n",    "n",    1000., -100., 5000.);
    RooRealVar G000("G000", "G000",  0.5);
    RooRealVar G001("G001", "G001",  g1, -5., 5.);
    RooRealVar G002("G002", "G002",  g2, -5., 5.);
    RooRealVar G003("G003", "G003", 0.0);//, -1., 1.);
    RooRealVar G004("G004", "G004", 0.0);//, -1., 1.);
    RooRealVar    a("a",    "a",    va);
    RooRealVar    b("b",    "b",    vb);
    RooRealVar    n("n",    "n",    1000., -100., 5000.);

    RooB2Kll pdf("pdf", "pdf", *cosTheta, G000, G001, G002, G003, G004, a, b, n);

    RooFitResult * fitresult = pdf.fitTo(*data,Save(kTRUE), Minos(kFALSE), NumCPU(4), SumW2Error(kTRUE));
    
    RooPlot * frame = cosTheta->frame();
    data->plotOn(frame);
    pdf.plotOn(frame);

    if(fitresult->minNll() == fitresult->minNll() && fitresult->minNll()>0 ) {
    fout << i << "\t" << a.getVal() << "\t" << b.getVal() << "\t" << fitresult->minNll() << "\t" << fitresult->status() << "\t" 
         << G000.getVal() << "\t" << G001.getVal() << "\t" << G002.getVal() << "\t" << G003.getVal() << "\t" << G004.getVal() << endl;
    }

    gROOT->ProcessLine(".x ~/lhcb/lhcbStyle.C");
    TCanvas c("fit","fit", 800, 800);
    frame->Draw();
    c.SaveAs("fits/fit"+TString(namestr)+".png");
    c.SaveAs("fits/fit"+TString(namestr)+".pdf");
}

///
/// Fills vector with floating pars names
///
void Utils::getParameters(const RooFitResult &result, std::vector<TString> &names){
	RooArgList pars           = result.floatParsFinal();
	TIterator * it            = pars.createIterator();
	while(RooRealVar* p = (RooRealVar*) it->Next()){
		names.push_back(TString(p->GetName()));
	}
};

void runFit(RooAbsPdf *pdf, RooDataSet *data, double *NLL, int *stat_t, int MaxTries, int mhLow, int mhHigh){

	int ntries=0;
	int stat=1;
	double minnll=10e8;
	while (stat!=0){
	  if (ntries>=MaxTries) break;
	  RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),Range(mhLow,mhHigh));
	  //RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),Range(85,110));
	  //RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),SumW2Error(kTRUE)
          stat = fitTest->status();
	  minnll = fitTest->minNll();
	  ntries++; 
	}
	*stat_t = stat;
	*NLL = minnll;
}

RooFitResult* FitHistWithCBShape(TH1* h, double mass)
{
  RooFitResult *res = 0;

  double avg = 0.97*mass;
  double rms = 0.05*mass;
  double xmin = 0.7*mass;
  // double xmax = 1.25*mass;
  // double xmin = 0.80*mass;
  double xmax = 1.25*mass;

  RooRealVar E("E","Normalised Resonance Mass", xmin, xmax, "GeV");
  RooDataHist dataSet("dataSet", "dataSet", E, h ); 

  RooPlot* frame = E.frame();
  dataSet.plotOn(frame);

  RooRealVar mean("mean","mean", avg, 0.92*mass, 1.2*mass) ;
  RooRealVar sigma("sigma","width", rms, 0.001*mass, 0.2*mass);
  RooRealVar alphaLow("alphaLow","alpha", 0.8, 0.0, 3.0);
  RooRealVar nLow("nLow","n", 20, 0.00001, 1000.0);
  RooRealVar alphaHigh("alphaHigh","alphaHigh", 2.0, 0.0, 3.0);
  RooRealVar nHigh("nHigh","nHigh", 100.0, 0.0000000001, 1000.0);

  RooFormulaVar minusE("minusE","Negate the data set","-1*E", E) ;
  RooFormulaVar minusMean("minusMean","Negate the mean","-1*mean", mean) ;


  RooCBShape CBallLow("CBallLow","Crystal Ball for lower tail", E, mean, sigma, alphaLow, nLow);
  RooCBShape CBallHigh("CBallHigh","Crystal Ball for highr tail", minusE, minusMean, sigma, alphaHigh, nHigh);  

  RooRealVar cb1frac("cb1frac","fraction of cb1", 0.5, 0, 1.0) ;

  RooAddPdf sum("sum","cb1+cb2",RooArgList(CBallLow, CBallHigh), cb1frac) ;

  res = sum.fitTo(dataSet, Save() );
  sum.plotOn(frame);
	
  frame->Draw();

  double sMean = ((RooAbsReal*) res->floatParsFinal().find("mean"))->getVal();
  double sSigma = ((RooAbsReal*) res->floatParsFinal().find("sigma"))->getVal();
  double sAlphaLow = ((RooAbsReal*) res->floatParsFinal().find("alphaLow"))->getVal();
  double sAlphaHigh = ((RooAbsReal*) res->floatParsFinal().find("alphaHigh"))->getVal();
  double sNLow = ((RooAbsReal*) res->floatParsFinal().find("nLow"))->getVal();
  double sNHigh = ((RooAbsReal*) res->floatParsFinal().find("nHigh"))->getVal();
  double sFrac = ((RooAbsReal*) res->floatParsFinal().find("cb1frac"))->getVal();

  Masses.push_back(mass); Means->push_back(sMean); Sigmas->push_back(sSigma);
  alphaHights->push_back(sAlphaHigh), alphaLows->push_back(sAlphaLow);
  nHighs->push_back(sNHigh), nLows->push_back(sNLow), fracs->push_back(sFrac);


  return res;

}

void JeffreysPriorDemo(){
  RooWorkspace w("w");
  w.factory("Uniform::u(x[0,1])");
  w.factory("mu[100,1,200]");
  w.factory("ExtendPdf::p(u,mu)");

  //  w.factory("Poisson::pois(n[0,inf],mu)");

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
  //  RooDataHist* asimov2 = w.pdf("pois")->generateBinned(*w.var("n"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  w.defineSet("poi","mu");
  w.defineSet("obs","x");
  //  w.defineSet("obs2","n");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",10);

  //  JeffreysPrior pi2("jeffreys2","jeffreys",*w.pdf("pois"),*w.set("poi"),*w.set("obs2"));

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1./sqrt(mu)",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  //  pi.plotOn(plot, Normalization(1,RooAbsReal::Raw),Precision(.1));
  pi.plotOn(plot);
  //  pi2.plotOn(plot,LineColor(kGreen),LineStyle(kDotted));
  test->plotOn(plot,LineColor(kRed));
  plot->Draw();

}

//_________________________________________________
void TestJeffreysGaussMean(){
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
  w.factory("n[10,.1,200]");
  w.factory("ExtendPdf::p(g,n)");
  w.var("sigma")->setConstant();
  w.var("n")->setConstant();

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  //  w.defineSet("poi","mu,sigma");
  w.defineSet("poi","mu");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}