BayesianAnalysisToolkit: BATCalculator.cxx Source File

Go to the documentation of this file.
00001 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Stefan A. Schmitz
00002 
00003 #include <TAxis.h>
00004 #include <TFile.h>
00005 #include <TCanvas.h>
00006 #include <TVector.h>
00007 
00008 #include <RooAbsFunc.h>
00009 //#include <RooAbsReal.h>
00010 #include <RooRealVar.h>
00011 //#include <RooArgSet.h>
00012 #include <RooBrentRootFinder.h>
00013 #include <RooFormulaVar.h>
00014 #include <RooGenericPdf.h>
00015 //#include <RooPlot.h>
00016 #include <RooProdPdf.h>
00017 #include <RooDataHist.h>
00018 //#include <RooAbsPdf.h>
00019 #include <RooHistPdf.h>
00020 
00021 #include <RooStats/MarkovChain.h>
00022 
00023 #include "../../BAT/BCLog.h"
00024 #include "../../BAT/BCAux.h"
00025 #include "../../BAT/BCH1D.h"
00026 
00027 #include "BCRooInterface.h"
00028 #include "BATCalculator.h"
00029 
00030 #include <algorithm>
00031 #include <cassert>
00032 
00033 
00034 ClassImp(RooStats::BATCalculator)
00035 
00036 namespace RooStats {
00037 
00038 // ---------------------------------------------------------
00039 BATCalculator::BATCalculator()
00040    : fData(0)
00041    , fPdf(0)
00042    , fPrior(0)
00043    , fProductPdf(0)
00044    , fLogLike(0)
00045    , fLikelihood(0)
00046    , fIntegratedLikelihood(0)
00047    , fPosteriorPdf(0)
00048    , fLower(0)
00049    , fUpper(0)
00050    , fBrfPrecision(0.00005)
00051    , fValidInterval(false)
00052    , fValidMCMCInterval(false)
00053    , fSize(0.05)
00054    , fLeftSideFraction(0.5)
00055 {
00056    // default constructor
00057    _myRooInterface = new BCRooInterface();
00058 }
00059 
00060 // ---------------------------------------------------------
00061 BATCalculator::BATCalculator( /* const char * name,  const char * title, */
00062                               RooAbsData & data,
00063                               RooAbsPdf  & pdf,
00064                               RooArgSet  & POI,
00065                               RooAbsPdf  & prior,
00066                               RooArgSet  * params,
00067                               bool fillChain)
00068    : fData(&data)
00069    , fPdf(&pdf)
00070    , fPOI(POI)
00071    , fPrior(&prior)
00072    , fparams(params)
00073    , fProductPdf(0)
00074    , fLogLike(0)
00075    , fLikelihood(0)
00076    , fIntegratedLikelihood(0)
00077    , fPosteriorPdf(0)
00078    , fLower(0)
00079    , fUpper(0)
00080    , fBrfPrecision(0.00005)
00081    , fValidInterval(false)
00082    , fValidMCMCInterval(false)
00083    , _nMCMC(1000000)
00084    , fSize(0.05)
00085    , fLeftSideFraction(0.5)
00086    //, TNamed( TString(name), TString(title) )
00087 {
00088    // constructor
00089    //if (nuisanceParameters)
00090    //   fNuisanceParameters.add(*nuisanceParameters);
00091    _myRooInterface = new BCRooInterface("BCRooInterfaceForBAT",fillChain);
00092 }
00093 
00094 // ---------------------------------------------------------
00095 BATCalculator::BATCalculator( RooAbsData & data, ModelConfig & model, bool fillChain)
00096    : fData(&data)
00097    , fPdf(model.GetPdf())
00098    , fPOI(*model.GetParametersOfInterest())
00099    , fPrior(model.GetPriorPdf())
00100    , fparams(model.GetNuisanceParameters())
00101    , fProductPdf(0)
00102    , fLogLike(0)
00103    , fLikelihood(0)
00104    , fIntegratedLikelihood(0)
00105    , fPosteriorPdf(0)
00106    , fLower(0)
00107    , fUpper(0)
00108    , fBrfPrecision(0.00005)
00109    , fValidInterval(false)
00110    , fValidMCMCInterval(false)
00111    , _nMCMC(1000000)
00112    , fSize(0.05)
00113    , fLeftSideFraction(0.5) 
00114 {
00115    std::cout << "BATCalculator calling constructor ..." << std::endl;
00116    // constructor from Model Config
00117    //  SetModel(model);
00118    _myRooInterface = new BCRooInterface("BCRooInterfaceForBAT",fillChain);
00119    std::cout << "BATCalculator constructed" << std::endl;
00120 }
00121 
00122 // ---------------------------------------------------------
00123 BATCalculator::~BATCalculator()
00124 {
00125    std::cout << "BATCalculator calling destructor ..." << std::endl;
00126    // destructor
00127    ClearAll();
00128    delete _myRooInterface;
00129 }
00130 
00131 // ---------------------------------------------------------
00132 void BATCalculator::ClearAll() const
00133 {
00134    // clear cached pdf objects
00135    if (fProductPdf)
00136       delete fProductPdf;
00137    if (fLogLike)
00138       delete fLogLike;
00139    if (fLikelihood)
00140       delete fLikelihood;
00141    if (fIntegratedLikelihood)
00142       delete fIntegratedLikelihood;
00143    if (fPosteriorPdf)
00144       delete fPosteriorPdf;
00145    fPosteriorPdf = 0;
00146    fProductPdf = 0;
00147    fLogLike = 0;
00148    fLikelihood = 0;
00149    fIntegratedLikelihood = 0;
00150    fLower = 0;
00151    fUpper = 0;
00152    fValidInterval = false;
00153    fValidMCMCInterval = false;
00154 }
00155 
00156 // ---------------------------------------------------------
00157 void BATCalculator::SetModel(const ModelConfig & /* model */)
00158 {
00159    // set the model
00160    //fPdf = model.GetPdf();
00161    //fPriorPOI =  model.GetPriorPdf();
00162    // assignment operator = does not do a real copy the sets (must use add method)
00163    //fPOI.removeAll();
00164    //*fparams.removeAll();
00165    //if (model.GetParametersOfInterest())
00166    //   fPOI.add( *(model.GetParametersOfInterest()) );
00167    //if (model.GetNuisanceParameters())
00168    //   fNuisanceParameters.add( *(model.GetNuisanceParameters() ) );
00169 
00170    // invalidate the cached pointers
00171    //ClearAll();
00172 }
00173 
00174 // ---------------------------------------------------------
00175 RooArgSet * BATCalculator::GetMode(RooArgSet * /* parameters */) const
00176 {
00179   // Should cover multi-dimensional cases...
00180   // How do we do if there are points that are equi-probable?
00181 
00182   return 0;
00183 }
00184 
00185 // ---------------------------------------------------------
00186 //returns posterior with respect to first entry in the POI ArgSet
00187 RooAbsPdf * BATCalculator::GetPosteriorPdf1D() const
00188 {
00189    const char * POIname = fPOI.first()->GetName();
00190    return GetPosteriorPdf1D(POIname);
00191 }
00192 
00193 // ---------------------------------------------------------
00194 // returns posterior with respect to provided name in the POI ArgSet
00195 RooAbsPdf * BATCalculator::GetPosteriorPdf1D(const char * POIname) const
00196 {
00197 
00198    // run some sanity checks
00199    if (!fPdf ) {
00200      std::cerr << "BATCalculator::GetPosteriorPdf - missing pdf model" << std::endl;
00201      return 0;
00202    }
00203 
00204    if (!fPrior) {
00205       std::cerr << "BATCalculator::GetPosteriorPdf - missing prior pdf" << std::endl;
00206    }
00207 
00208    if (fPOI.getSize() == 0) {
00209      std::cerr << "BATCalculator::GetPosteriorPdf - missing parameter of interest" << std::endl;
00210      return 0;
00211    }
00212 
00213    if (fPOI.getSize() > 1) {
00214       std::cerr << "BATCalculator::GetPosteriorPdf - current implementation works only on 1D intervals" << std::endl;
00215       return 0;
00216    }
00217 
00218     //initialize RooInterface object
00219    _myRooInterface->Initialize(*fData,*fPdf,*fPrior,fparams,fPOI);
00220    _myRooInterface->MCMCSetNIterationsRun(_nMCMC);
00221    _myRooInterface->MarginalizeAll();
00222    _myRooInterface->FindMode();
00223    BCParameter * myPOI = _myRooInterface->GetParameter(POIname);
00224    BCH1D * myPosterior =_myRooInterface->GetMarginalized(myPOI);
00225    TH1D * posteriorTH1D = myPosterior->GetHistogram();
00226    _posteriorTH1D = static_cast<TH1D *>(posteriorTH1D->Clone("_posteriorTH1D"));
00227    RooDataHist * posteriorRooDataHist = new RooDataHist("posteriorhist","", fPOI,posteriorTH1D);
00228    fPosteriorPdf = new RooHistPdf("posteriorPdf","",fPOI,*posteriorRooDataHist);
00229 
00230    /*
00231    // plots for debugging
00232    TFile * debugfile = new TFile( "debug_posterior.root" ,"RECREATE");
00233    if ( debugfile->IsOpen() )
00234       std::cout << "File debug_posterior opened successfully" << std::endl;
00235    debugfile->cd();
00236    //posteriorTH1D->Write("posteriorTH1D");
00237    //posteriorRooDataHist->Write("posteriorRooDataHist");
00238    fPosteriorPdf->Write("fPosteriorPdf");
00239    TCanvas myCanvas("canvasforfposterior","canvas for fposterior");
00240    myCanvas.cd();
00241    //fPosteriorPdf->Draw();
00242    //myCanvas.Write();
00243    RooRealVar nSig("nSig","",0.0001,200.);
00244    RooPlot * pl = nSig.frame();
00245    fPosteriorPdf->plotOn(pl);
00246    pl->Draw();
00247    pl->Write();
00248    myCanvas.Write();
00249    debugfile->Close();
00250 
00251    // just for testing: create plots with BAT
00252    char * outputFile = "bat_plots_debugging.ps";
00253    _myRooInterface -> PrintAllMarginalized(outputFile);
00254    //*/
00255 
00256    return fPosteriorPdf; // is of type RooAbsPdf*
00257 }
00258 
00259 // ---------------------------------------------------------
00260 // return a RooPlot with the posterior PDF and the credibility region
00261 RooPlot * BATCalculator::GetPosteriorPlot1D() const
00262 {
00263 
00264    if (!fPosteriorPdf){
00265       std::cout << "BATCalculator::GetPosteriorPlot1D:"
00266                 << "Warning : posterior not available" << std::endl;
00267       GetPosteriorPdf1D();
00268    }
00269    if ((!fValidInterval) && (!fValidMCMCInterval)){
00270       std::cout << "BATCalculator::GetPosteriorPlot1D:"
00271                 << "Warning : interval not available" << std::endl;
00272       GetInterval1D();
00273    }
00274 
00275    RooAbsRealLValue * poi = dynamic_cast<RooAbsRealLValue *>( fPOI.first() );
00276    assert(poi);
00277 
00278    RooPlot* plot = poi->frame();
00279 
00280    plot->SetTitle(TString("Posterior probability of parameter \"")+TString(poi->GetName())+TString("\""));
00281    fPosteriorPdf->plotOn(plot,RooFit::Range(fLower,fUpper,kFALSE),RooFit::VLines(),RooFit::DrawOption("F"),RooFit::MoveToBack(),RooFit::FillColor(kGray));
00282    fPosteriorPdf->plotOn(plot);
00283    plot->GetYaxis()->SetTitle("posterior probability");
00284 
00285    return plot;
00286 }
00287 
00288 // ---------------------------------------------------------
00289 // returns central interval for first POI in the POI ArgSet
00290 SimpleInterval * BATCalculator::GetInterval1D() const
00291 {
00292    const char * POIname = fPOI.first()->GetName();
00293    return GetInterval1D(POIname); //is of type RooAbsPdf *
00294 }
00295 
00296 
00297 // ---------------------------------------------------------
00298 // returns central interval for the defined POI in the POI ArgSet -> test code because orginal version is not working anymore
00299 SimpleInterval * BATCalculator::GetInterval1D(const char * POIname) const
00300 {
00301    //const char * POIname = fPOI.first()->GetName();
00302 
00303    if (fValidInterval)
00304    std::cout << "BATCalculator::GetShortestInterval1D:"
00305              << "Warning : recomputing interval for the same CL and same model" << std::endl;
00306 
00307    // get pointer to selected parameter of interest
00308    RooRealVar * poi = dynamic_cast<RooRealVar*>( fPOI.find(POIname) );
00309    assert(poi);
00310 
00311    // get pointer to poserior pdf
00312    if (!fPosteriorPdf)
00313       fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf1D();
00314    //RooAbsReal * cdf = fPosteriorPdf->createCdf(fPOI,RooFit::ScanParameters(100,2));
00315 
00316    // range of poi
00317    Double_t minpoi = poi->getMin();
00318    Double_t maxpoi = poi->getMax();
00319 
00320    // bin number of histogram for posterior
00321    Int_t stepnumber = _posteriorTH1D->GetNbinsX();
00322    std::cout << "stepnumber is: " << stepnumber << std::endl;
00323 
00324    // width of one bin in histogram of posterior
00325    Double_t stepsize = (maxpoi-minpoi)/stepnumber;
00326    std::cout << "stepsize is: " << stepsize << std::endl;
00327 
00328    // pair: first entry: bin number , second entry: value of posterior
00329    std::vector< std::pair< Int_t,Double_t > > posteriorVector;
00330 
00331    // for normalization
00332    Double_t histoIntegral = 0;
00333    // give posteriorVector the right length
00334    posteriorVector.resize(stepnumber);
00335 
00336    // see in BayesianCalculator for details about this "feature"
00337    // FB: comment in because 'unused variable'
00338    //   Double_t tmpVal = poi->getVal();
00339 
00340   // initialize elements of posteriorVector
00341    int i = 0;
00342    std::vector< std::pair< Int_t,Double_t > >::iterator vecit = posteriorVector.begin();
00343    std::vector< std::pair< Int_t,Double_t > >::iterator vecit_end = posteriorVector.end();
00344    for( ; vecit != vecit_end ; ++vecit) {
00345       poi->setVal(poi->getMin()+i*stepsize);
00346       posteriorVector[i] = std::make_pair(i, _posteriorTH1D->GetBinContent(i+1) ); // hope this is working, +1 necessary, because GetBinContent(0) returns the underflow bin
00347       histoIntegral+=_posteriorTH1D->GetBinContent(i); // better to get this directly from the histogram ?!
00348       //std::cout << "pair with binnumber: " << i << " and postriorprob: " << _posteriorTH1D->GetBinContent(i+1) << std::endl;
00349       i++;
00350    }
00351 
00352    std::cout << "histoIntegral is: " << histoIntegral << std::endl;
00353 
00354    double lowerProbLim = (1.-ConfidenceLevel())*fLeftSideFraction;
00355    double upperProbLim = 1. -( (1.-ConfidenceLevel())*(1.-fLeftSideFraction) );
00356    std::cout << "lowerProbLim is: " << lowerProbLim  << "upperProbLim is: " << histoIntegral << std::endl;
00357 
00358    fLower = -1.;
00359    fUpper = -1.;
00360    bool lowerlimset = false;
00361    bool upperlimset = false;
00362 
00363 
00364    if(fLeftSideFraction == 0){
00365       fLower = minpoi;
00366       lowerlimset = true;
00367    }
00368 
00369    if(fLeftSideFraction == 1){
00370       fUpper = maxpoi;
00371       upperlimset = true;
00372    }
00373 
00374    // keep track of integrated posterior in the interval
00375    Double_t integratedposterior = 0.;
00376 
00377    i = 0;
00378    vecit = posteriorVector.begin();
00379    for( ; vecit != vecit_end ; ++vecit) {
00380       integratedposterior += posteriorVector[i].second;
00381       if( (lowerlimset != true) && ((integratedposterior/histoIntegral) >= lowerProbLim) ){
00382          fLower = poi->getMin()+i*stepsize;
00383          lowerlimset = true;
00384       }
00385       if( (upperlimset != true) && ((integratedposterior/histoIntegral) >= upperProbLim) ){
00386          fUpper = poi->getMin()+i*stepsize;
00387          upperlimset = true;
00388          break;
00389       }
00390       i++;
00391    }
00392 
00393    TString interval_name = TString("CentralBayesianInterval_a") + TString(this->GetName());
00394    SimpleInterval * interval = new SimpleInterval(interval_name,*poi,fLower,fUpper,ConfidenceLevel());
00395    interval->SetTitle("SimpleInterval from BATCalculator");
00396 
00397    return interval;
00398 
00399 }
00400 
00401 // ---------------------------------------------------------
00402 // returns central interval for requested POI -> no longer working since root v5_32
00403 // SimpleInterval * BATCalculator::GetInterval1Dv0(const char * POIname) const
00404 // {
00405 //   /// returns a SimpleInterval with lower and upper bounds on the
00406 //   /// parameter of interest. Applies the central ordering rule to
00407 //   /// compute the credibility interval. Covers only the case with one
00408 //   /// single parameter of interest
00409 // 
00410 //    if (fValidInterval)
00411 //       std::cout << "BATCalculator::GetInterval1D:"
00412 //                 << "Warning : recomputing interval for the same CL and same model" << std::endl;
00413 // 
00414 //    RooRealVar * poi = dynamic_cast<RooRealVar *>( fPOI.find(POIname) );
00415 //    assert(poi);
00416 // 
00417 //    if (!fPosteriorPdf)
00418 //       fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf1D();
00419 // 
00420 //    Double_t* myarr = new Double_t[1];
00421 //    myarr[0]= 0.001;
00422 //    std::cout << "fPosteriorPdf: " << fPosteriorPdf->maxVal(1) << std::endl;
00423 //    std::cout << "fPosteriorPdf: " << fPosteriorPdf->getVal() << std::endl;
00424 // //   std::cout << "fPosteriorPdf: " << fPosteriorPdf->evaluate() << std::endl;
00425 // 
00426 //    RooAbsReal * cdf = fPosteriorPdf->createCdf(fPOI,RooFit::ScanParameters(100,2));
00427 //    //RooAbsReal* cdf = fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf());
00428 // 
00429 //    RooAbsFunc * cdf_bind = cdf->bindVars(fPOI,&fPOI);
00430 //    std::cout << "cdf_bind: " << std::endl;
00431 //    std::cout << "cdf_bind: " << cdf_bind->isValid() << std::endl;
00432 //    std::cout << "cdf_bind: " << cdf_bind->getDimension() << std::endl;
00433 //    //Double_t test = (*_function)(&a) - value;
00434 //    std::cout << "cdf_bind: " << (*cdf_bind)(myarr) << std::endl;
00435 //    std::cout << "cdf_bind: " << cdf_bind->getMaxLimit(1) << std::endl;
00436 // 
00437 //    RooBrentRootFinder brf(*cdf_bind);
00438 //    brf.setTol(fBrfPrecision); // set the brf precision
00439 // 
00440 //    double tmpVal = poi->getVal();  // patch used because findRoot changes the value of poi
00441 // 
00442 //    double y = fSize*fLeftSideFraction; //adjust lower tail prob. according to fLeftSideFraction
00443 // 
00444 //    brf.findRoot(fLower,poi->getMin(),poi->getMax(),y);
00445 // 
00446 //    y = 1.-(fSize*(1.-fLeftSideFraction) ); //adjust upper tail prob. according to fLeftSideFraction
00447 // 
00448 //    bool ret = brf.findRoot(fUpper,poi->getMin(),poi->getMax(),y);
00449 //    std::cout << " brf.findRoot: " << fUpper << " " << poi->getMin() << " " << poi->getMax() << " " << y << std::endl;
00450 //    if (!ret)
00451 //       std::cout << "BATCalculator::GetInterval1D: Warning:"
00452 //                 << "Error returned from Root finder, estimated interval is not fully correct"
00453 //                 << std::endl;
00454 // 
00455 //    poi->setVal(tmpVal); // patch: restore the original value of poi
00456 // 
00457 //    delete cdf_bind;
00458 //    delete cdf;
00459 //    fValidInterval = true;
00460 //    fConnectedInterval = true;
00461 // 
00462 //    TString interval_name = TString("CentralBayesianInterval_a") + TString(this->GetName());
00463 //    SimpleInterval * interval = new SimpleInterval(interval_name,*poi,fLower,fUpper,ConfidenceLevel());
00464 //    interval->SetTitle("SimpleInterval from BATCalculator");
00465 // 
00466 //    return interval;
00467 // }
00468 
00469 // ---------------------------------------------------------
00470 SimpleInterval * BATCalculator::GetShortestInterval1D() const
00471 {
00472    const char * POIname = fPOI.first()->GetName();
00473    bool checkConnected = true;
00474    return GetShortestInterval1D(POIname, checkConnected);
00475 }
00476 
00477 // ---------------------------------------------------------
00478 // returns a SimpleInterval with lower and upper bounds on the
00479 // parameter of interest. Applies the shortest interval rule to
00480 // compute the credibility interval. The resulting interval is not necessarily connected.
00481 // Methods for constructing the shortest region with respect given CL are now available in
00482 // different places (here; MCMCCalculator, BAT)-> this might require cleaning at some point.
00483 // Result is approximate as CL is not reached exactly (due to finite bin number/bin size in posterior histogram)-> make CL/interval a bit smaller or bigger ?
00484 SimpleInterval * BATCalculator::GetShortestInterval1D(const char * POIname, bool & checkConnected) const
00485 {
00486 
00487    if (fValidInterval)
00488       std::cout << "BATCalculator::GetShortestInterval1D:"
00489                 << "Warning : recomputing interval for the same CL and same model" << std::endl;
00490 
00491    // get pointer to selected parameter of interest
00492    RooRealVar * poi = dynamic_cast<RooRealVar*>( fPOI.find(POIname) );
00493    assert(poi);
00494 
00495    // get pointer to poserior pdf
00496    if (!fPosteriorPdf)
00497       fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf1D();
00498    //RooAbsReal * cdf = fPosteriorPdf->createCdf(fPOI,RooFit::ScanParameters(100,2));
00499 
00500    // range of poi
00501    Double_t minpoi = poi->getMin();
00502    Double_t maxpoi = poi->getMax();
00503 
00504    // bin number of histogram for posterior
00505    Int_t stepnumber = _posteriorTH1D->GetNbinsX();
00506    std::cout << "stepnumber is: " << stepnumber << std::endl;
00507 
00508    // width of one bin in histogram of posterior
00509    Double_t stepsize = (maxpoi-minpoi)/stepnumber;
00510    std::cout << "stepsize is: " << stepsize << std::endl;
00511 
00512    // pair: first entry: bin number , second entry: value of posterior
00513    std::vector< std::pair< Int_t,Double_t > > posteriorVector;
00514 
00515    // for normalization
00516    Double_t histoIntegral = 0;
00517    // give posteriorVector the right length
00518    posteriorVector.resize(stepnumber);
00519 
00520    // see in BayesianCalculator for details about this "feature"
00521    Double_t tmpVal = poi->getVal();
00522 
00523    // initialize elements of posteriorVector
00524    int i = 0;
00525    std::vector< std::pair< Int_t,Double_t > >::iterator vecit = posteriorVector.begin();
00526    std::vector< std::pair< Int_t,Double_t > >::iterator vecit_end = posteriorVector.end();
00527    for( ; vecit != vecit_end ; ++vecit) {
00528       poi->setVal(poi->getMin()+i*stepsize);
00529       posteriorVector[i] = std::make_pair(i, _posteriorTH1D->GetBinContent(i+1) ); // hope this is working, +1 necessary, because GetBinContent(0) returns the underflow bin
00530       histoIntegral+=_posteriorTH1D->GetBinContent(i); // better to get this directly from the histogram ?!
00531       //std::cout << "pair with binnumber: " << i << " and postriorprob: " << _posteriorTH1D->GetBinContent(i+1) << std::endl;
00532       i++;
00533    }
00534 
00535    std::cout << "histoIntegral is: " << histoIntegral << std::endl;
00536 
00537    // sort posteriorVector with respect to posterior pdf
00538    std::sort(posteriorVector.begin(), posteriorVector.end(), sortbyposterior);
00539 
00540    // keep track of integrated posterior in the interval
00541    Double_t integratedposterior = 0.;
00542 
00543    // keep track of lowerLim and upperLim
00544    Double_t lowerLim=posteriorVector.size();
00545    Double_t upperLim=0;
00546 
00547    // store the bins in the intervall
00548    std::vector<bool> inInterval;
00549    inInterval.resize(posteriorVector.size());
00550 
00551    // set all values in inInterval to false
00552    for (unsigned int k = 0; k < inInterval.size(); k++)
00553       inInterval[k] = false;
00554 
00555    unsigned int j = 0;
00556    // add bins to interval while CL not reached
00557    //std::cout << "integratedposterior: " << integratedposterior << "histoIntegral: " << histoIntegral << std::endl;
00558    //std::cout << " integratedposterior/histoIntegral : " << integratedposterior/histoIntegral << std::endl;
00559    //std::cout << " 1-fSize : " << 1-fSize << std::endl;
00560    //std::cout << " posteriorVector.size(): " << posteriorVector.size() << std::endl;
00561    while(((integratedposterior/histoIntegral) < (1-fSize)) && (j < posteriorVector.size())) {
00562       //std::cout << "j is: " << j << " , int prob: " << integratedposterior/histoIntegral << std::endl;
00563       integratedposterior+=posteriorVector[j].second;
00564 
00565       //std::cout << "bin number: " << posteriorVector[j].first << " with posterior prob.: " << posteriorVector[j].second << std::endl;
00566 
00567       // update vector with bins included in the interval
00568       inInterval[posteriorVector[j].first] = true;
00569 
00570       if(posteriorVector[j].first < lowerLim) {
00571          lowerLim = posteriorVector[j].first; // update lowerLim
00572          std::cout << "updating lower lim to: " << lowerLim << std::endl;
00573       }
00574       if(posteriorVector[j].first > upperLim) {
00575          upperLim = posteriorVector[j].first; // update upperLim
00576          std::cout << "updating upper lim to: " << upperLim << std::endl;
00577       }
00578 
00579       fLower = lowerLim * stepsize; // 
00580       fUpper = upperLim * stepsize;
00581       j++;
00582    }
00583 
00584    // initialize vector with interval borders
00585 
00586    bool runInside = false;
00587    for (unsigned int l = 0; l < inInterval.size(); l++) {
00588       if ( (runInside == false) && (inInterval[l] == true) ) {
00589          _intervalBorders1D.push_back(static_cast<double>(l)* stepsize);
00590          runInside = true;
00591       }
00592       if ( ( runInside == true) && (l < (inInterval.size()-1) ) && (inInterval[l+1] == false) ) {
00593          _intervalBorders1D.push_back(static_cast<double>(l)* stepsize);
00594          runInside = false;
00595       }
00596       if ( ( runInside == true) && (l == (inInterval.size()-1)) ) {
00597          _intervalBorders1D.push_back(static_cast<double>(l)* stepsize);
00598       }
00599    }
00600 
00601    // check if the intervall is connected
00602    if(checkConnected) {
00603       if (_intervalBorders1D.size() > 2) {
00604          fConnectedInterval = false;
00605       }
00606       else {
00607          fConnectedInterval = true;
00608       }
00609    }
00610 
00611    poi->setVal(tmpVal); // patch: restore the original value of poi
00612 
00613    fValidInterval = true;
00614 
00615    TString interval_name = TString("ShortestBayesianInterval_a") + TString(this->GetName());
00616    SimpleInterval * interval = new SimpleInterval(interval_name,*poi,fLower,fUpper,ConfidenceLevel());
00617    interval->SetTitle("Shortest SimpleInterval from BATCalculator");
00618 
00619    //if(assert)
00620    return interval;
00621 }
00622 
00623 // ---------------------------------------------------------
00624 MCMCInterval* BATCalculator::GetInterval() const{
00625 
00626    // run some sanity checks
00627    if (!fPdf ) {
00628      std::cerr << "BATCalculator::GetInterval - missing pdf model" << std::endl;
00629      return 0;
00630    }
00631 
00632    if (!fPrior) {
00633       std::cerr << "BATCalculator::GetInterval - missing prior pdf" << std::endl;
00634    }
00635 
00636    if (fPOI.getSize() == 0) {
00637      std::cerr << "BATCalculator::GetInterval - missing parameter of interest" << std::endl;
00638      return 0;
00639    }
00640 
00641    if (!fPosteriorPdf){
00642       //initialize RooInterface object
00643       _myRooInterface->Initialize(*fData,*fPdf,*fPrior,fparams,fPOI);
00644       _myRooInterface->MCMCSetNIterationsRun(_nMCMC);
00645       _myRooInterface->MarginalizeAll();
00646       _myRooInterface->FindMode();
00647    }
00648 
00649    MarkovChain * roostats_chain = GetBCRooInterface()->GetRooStatsMarkovChain();
00650    MCMCInterval * mcmcInterval = new MCMCInterval("roostatsmcmcinterval", fPOI , *roostats_chain);
00651    //MCMCInterval * mcmcInterval = new MCMCInterval("roostatsmcmcinterval", *(GetBCRooInterface()->GetArgSetForMarkovChain()) , *roostats_chain);
00652    mcmcInterval->SetUseKeys(false);
00653    mcmcInterval->SetConfidenceLevel(1.-fSize);
00654    fValidMCMCInterval = true;
00655    return mcmcInterval;
00656 }
00657 
00658 // ---------------------------------------------------------
00659 void BATCalculator::SetNumBins(const char * parname, int nbins)
00660 {
00661    _myRooInterface->SetNumBins(parname, nbins);
00662 }
00663 
00664 void BATCalculator::SetNumBins(int nbins)
00665 {
00666    _myRooInterface->SetNumBins(nbins);
00667 }
00668 
00669 void BATCalculator::SetLeftSideTailFraction(Double_t leftSideFraction ){
00670    if( (leftSideFraction >= 0.) && (leftSideFraction <= 1.) ){
00671       fLeftSideFraction = leftSideFraction;
00672    }
00673    else{
00674       std::cout << "BATCalculator::SetLeftSideTailFraction(Double_t leftSideFraction ) - value needs to be in the interval [0.,1.] to be meaningful, your value: " << leftSideFraction <<  " ,left side tail fraction has not been changed!" << std::endl;
00675    }
00676 
00677 } 
00678 
00679 
00680 // ---------------------------------------------------------
00681 Double_t BATCalculator::GetOneSidedUperLim()
00682 {
00683 //   double safeVal = fSize;
00684 //   fSize = safeVal/2.;
00685    std::cout << "calculating " << (1.-fSize/2) << "upper limit" << std::endl;
00686    return GetInterval1D()->UpperLimit();
00687 }
00688 
00689 /*
00690 // ---------------------------------------------------------
00691 int BATCalculator::GetNbins(const char * parname)
00692 {
00693    ;
00694 }
00695 */
00696 
00697 // ---------------------------------------------------------
00698 bool sortbyposterior(std::pair< Int_t,Double_t > pair1, std::pair< Int_t,Double_t > pair2)
00699 {
00700    return (pair1.second > pair2.second);
00701 }
00702 
00703 
00704 
00705 } // end namespace RooStats