• Main Page
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

BCRooInterface.cxx

Go to the documentation of this file.
00001 #include <RooGlobalFunc.h>
00002 #include <RooMsgService.h>
00003 #include <RooProdPdf.h>
00004 #include <RooWorkspace.h>
00005 #include <TFile.h>
00006 
00007 #include "../../BAT/BCMath.h"
00008 #include <iostream>
00009 
00010 #include "BCRooInterface.h"
00011 
00012 #include <RooUniform.h>
00013 
00014 //for testing
00015 #include "RooRealVar.h"
00016 #include "RooAbsReal.h"
00017 #include "RooRandom.h"
00018 
00019 
00020 
00021 
00022 // ---------------------------------------------------------
00023 void BCRooInterface::Initialize( RooAbsData& data,
00024              RooAbsPdf& model,
00025              RooAbsPdf& prior_trans,
00026              const RooArgSet* params,
00027              const RooArgSet& listPOI )
00028 {
00029 
00030   fData = &data;
00031   fModel = &model;
00032 
00033   // make the product of both priors to get the full prior probability function
00034   RooAbsPdf* prior_total = &prior_trans;
00035   if (prior_total!=0 ) {
00036     fPrior = prior_total;
00037   }
00038   else {
00039          std::cout << "No prior PDF: without taking action the program would crash\n";
00040          std::cout << "No prior PDF: adding dummy uniform pdf on the interval [0..1]\n";
00041          priorhelpvar = new RooRealVar("priorhelpvar","",0.0, 1.0 );
00042          _addeddummyprior = true;
00043          RooUniform* priorhelpdist = new RooUniform("priorhelpdist","", *priorhelpvar);
00044          fPrior = priorhelpdist;
00045   }
00046 
00047   std::cout << "Imported parameters:\n";
00048   fParams  = new RooArgList(listPOI);
00049   const RooArgSet* paramsTmp = params;
00050   if (paramsTmp!=0)
00051     fParams->add(*paramsTmp);
00052   fParams->Print("v");
00053 
00054   fParamsPOI  = new RooArgList(listPOI);
00055   // create the log-likelihood function
00056 //   fNll = new RooNLLVar("fNll","",*fModel,*fData,true/*extended*/);
00057 
00058   RooArgSet* constrainedParams = fModel->getParameters(*fData);
00059   fNll = fModel->createNLL(*fData, RooFit::Constrain(*constrainedParams) );
00060 
00061   DefineParameters();
00062 
00063   if(_fillChain){
00064       SetupRooStatsMarkovChain();
00065   }
00066 }
00067 
00068 // ---------------------------------------------------------
00069 void BCRooInterface::Initialize( const char* rootFile,
00070              const char* wsName,
00071              const char* dataName,
00072              const char* modelName,
00073              const char* priorName,
00074              const char* priorNuisanceName,
00075              const char* paramsName,
00076              const char* listPOIName )
00077 {
00078   // retrieve the RooFit inputs from the ROOT file
00079 
00080   /*
00081   // hard coded names in the workspace
00082   char* rootFile = "bat_workspace.root";
00083   char* wsName= "batWS";
00084   char* dataName= "data";
00085   char* modelName= "model";
00086   char* priorName= "priorPOI";
00087   char* priorNuisanceName= "priorNuisance";
00088   char* paramsName= "parameters";
00089   char* listPOIName= "POI";
00090   */
00091 
00092    std::cout << "Opening " << rootFile << std::endl;
00093    TFile* file = new TFile(rootFile);
00094    std::cout << "content :\n";
00095    file->ls();
00096 
00097    RooWorkspace* bat_ws = (RooWorkspace*) file->Get(wsName);
00098    bat_ws->Print("v");
00099 
00100    fData = (RooAbsData*) bat_ws->data(dataName);
00101    fModel = (RooAbsPdf*) bat_ws->function(modelName);
00102 
00103    // make the product of both priors to get the full prior probability function
00104    RooAbsPdf* priorPOI = (RooAbsPdf*) bat_ws->function(priorName);
00105    RooAbsPdf* priorNuisance = (RooAbsPdf*) bat_ws->pdf(priorNuisanceName);
00106    if (priorNuisance!=0 && priorPOI!=0) {
00107       fPrior = new RooProdPdf("fPrior","complete prior",*priorPOI,*priorNuisance);
00108    }
00109    else {
00110       if ( priorNuisance!=0 )
00111          fPrior=priorNuisance;
00112       else if ( priorPOI!=0 )
00113          fPrior = priorPOI;
00114       else
00115          std::cout << "No prior PDF: without taking action the program would crash\n";
00116          std::cout << "No prior PDF: adding dummy uniform pdf on the interval [0..1]\n";
00117          priorhelpvar = new RooRealVar("priorhelpvar","",0.0, 1.0 );
00118          _addeddummyprior = true;
00119          RooUniform* priorhelpdist = new RooUniform("priorhelpdist","", *priorhelpvar);
00120          fPrior = priorhelpdist;
00121    }
00122 
00123    std::cout << "Imported parameters:\n";
00124    fParams  = new RooArgList(*(bat_ws->set(listPOIName)));
00125    RooArgSet* paramsTmp = (RooArgSet*) bat_ws->set(paramsName);
00126    if (paramsTmp!=0) {
00127       fParams->add(*paramsTmp);
00128    }
00129    if (_addeddummyprior == true ) {
00130       fParams->add(*priorhelpvar);
00131    }
00132    fParams->Print("v");
00133 
00134    // create the log-likelihood function
00135    //fNll = new RooNLLVar("fNll","",*fModel,*fData,true/*extended*/);
00136    RooArgSet* constrainedParams = fModel->getParameters(*fData);
00137    fNll = fModel->createNLL(*fData, RooFit::Constrain(*constrainedParams) );
00138 
00139    file->Close();
00140 
00141    DefineParameters();
00142 }
00143 
00144 
00145 // ---------------------------------------------------------
00146 BCRooInterface::BCRooInterface() : BCModel()
00147 {   // default constructor
00148    _default_nbins = 500;
00149    _fillChain = false;
00150    _addeddummyprior = false;
00151 }
00152 
00153 // ---------------------------------------------------------
00154 BCRooInterface::BCRooInterface(const char* name, bool fillChain) : BCModel(name)
00155 {   // another constructor
00156    _default_nbins = 500;
00157    _fillChain = fillChain;
00158 }
00159 
00160 // ---------------------------------------------------------
00161 BCRooInterface::~BCRooInterface()
00162 {   // default destructor
00163    if(_fillChain){
00164       delete _roostatsMarkovChain;
00165    }
00166   //test begin
00167   //cout << "fIterationInterfacecount: " << fIterationInterfacecount << endl;
00168   //test end
00169 }
00170 
00171 // ---------------------------------------------------------
00172 void BCRooInterface::DefineParameters()
00173 {   // define for BAT the list of parameters, range and plot binning
00174 
00175   int nParams = fParams->getSize();
00176   for (int iParam=0; iParam<nParams; iParam++) {
00177     RooRealVar* ipar = (RooRealVar*) fParams->at(iParam);
00178     this->AddParameter(ipar->GetName(),ipar->getMin(),ipar->getMax());
00179     this->SetNbins(ipar->GetName(),_default_nbins);
00180     std::cout << "added parameter: " << ipar->GetName() << " defined in range [ " << ipar->getMin() << " - " << ipar->getMax() << " ] with number of bins: " << _default_nbins << " \n";
00181   }
00182 
00183    for(list< pair<const char*,int> >::iterator listiter = _nbins_list.begin(); listiter != _nbins_list.end(); listiter++){
00184       this->SetNbins((*listiter).first,(*listiter).second);
00185       std::cout << "adjusted parameter: " << (*listiter).first << " to number of bins: " << (*listiter).second << " \n";
00186    }
00187 
00188 }
00189 
00190 // ---------------------------------------------------------
00191 double BCRooInterface::LogLikelihood(const std::vector<double> & parameters)
00192 {   // this methods returns the logarithm of the conditional probability: p(data|parameters)
00193 
00194    // retrieve the values of the parameters to be tested
00195   int nParams = fParams->getSize();
00196   for (int iParam=0; iParam<nParams; iParam++) {
00197     RooRealVar* ipar = (RooRealVar*) fParams->at(iParam);
00198     ipar->setVal(parameters.at(iParam));
00199   }
00200 
00201   // compute the log of the likelihood function
00202   double logprob = -fNll->getVal();
00203   return logprob;
00204 }
00205 
00206 // ---------------------------------------------------------
00207 double BCRooInterface::LogAPrioriProbability(const std::vector<double> & parameters)
00208 {   // this method returs the logarithm of the prior probability for the parameters: p(parameters).
00209    // retrieve the values of the parameters to be tested
00210   int nParams = fParams->getSize();
00211   for (int iParam=0; iParam<nParams; iParam++) {
00212     RooRealVar* ipar = (RooRealVar*) fParams->at(iParam);
00213     ipar->setVal(parameters.at(iParam));
00214   }
00215   // compute the log of the prior function
00216   RooArgSet* tmpArgSet = new RooArgSet(*fParams);
00217   double prob = fPrior->getVal(tmpArgSet);
00218   delete tmpArgSet;
00219   if (prob<1e-300)
00220     prob = 1e-300;
00221   return log(prob);
00222 }
00223 
00224 // ---------------------------------------------------------
00225 void BCRooInterface::SetNumBins(const char * parname, int nbins)
00226 {
00227    for(list< pair<const char*,int> >::iterator listiter = _nbins_list.begin(); listiter != _nbins_list.end(); listiter++){
00228       if(!strcmp((*listiter).first, parname)){
00229          (*listiter).second = nbins;
00230          return;
00231       }
00232    }
00233    _nbins_list.push_back( make_pair(parname,nbins) );
00234 }
00235 
00236 // ---------------------------------------------------------
00237 void BCRooInterface::SetNumBins(int nbins)
00238 {
00239    _default_nbins = nbins;
00240 }
00241 
00242 // ---------------------------------------------------------
00243 void BCRooInterface::SetupRooStatsMarkovChain()
00244 {
00245    //RooArgSet * tempRooArgSetPointer = new RooArgSet(* fParams, "paramsMarkovChainPointer");
00246    //_parametersForMarkovChain = RooArgSet(* fParams, "paramsMarkovChainPointer");
00247    //fParams->snapshot(_parametersForMarkovChain);
00248    _parametersForMarkovChainPrevious.add(*fParamsPOI);
00249    _parametersForMarkovChainCurrent.add(*fParamsPOI);
00250 
00251    cout << "size of _parametersForMarkovChain: " << _parametersForMarkovChainCurrent.getSize() << endl;
00252    cout << "size of fParamsPOI: " << fParamsPOI->getSize() << endl;
00253    //cout << "size of tempRooArgSetPointer: " << tempRooArgSetPointer->getSize() << endl;
00254 
00255     _roostatsMarkovChain = new RooStats::MarkovChain();
00256    //test stuff begin
00257    //the following way of creating the MArkovChain object does not work!:
00258    //_roostatsMarkovChain = new RooStats::MarkovChain(_parametersForMarkovChain);
00259    //test stuff end
00260    cout << "setting up parameters for RooStats markov chain" << endl;
00261    _parametersForMarkovChainPrevious.writeToStream(cout, false);
00262 
00263    //initialize fPreviousStep, fCurrentStep, fVecWeights
00264    int nchains = MCMCGetNChains();
00265    for(int countChains = 1; countChains<=nchains ; countChains++ ) {
00266       double tempweight = 1.0;
00267       fVecWeights.push_back(tempweight);
00268       vector<double> tempvec;
00269       TIterator* setiter = fParamsPOI->createIterator();
00270       double tempval = 0.;
00271       while(setiter->Next()) {
00272          tempvec.push_back(tempval);
00273       }
00274       fPreviousStep.push_back(tempvec);
00275       fCurrentStep.push_back(tempvec);
00276    }
00277 
00278    fFirstComparison = true;
00279 
00280    //test stuff begin
00281    //var1 = new RooRealVar("var1","var1", 10., 0., 100.);
00282    //fParamsTest  = new RooArgList();
00283    //fParamsTest->add(*var1);
00284    //_parametersForMarkovChain_test.add(*fParamsTest);
00285    //fIterationInterfacecount = 0;
00286    //test stuff end
00287 }
00288 
00289 // ---------------------------------------------------------
00290 // to be added: add elements with higher weights if elements repeat to avoid unneccessarily long chains
00291 void BCRooInterface::MCMCIterationInterface()
00292 {
00293   //fIterationInterfacecount+=1;
00294 
00295   if(_fillChain){
00296       //cout << "Hei-ho running with fillChain!" << endl;
00297       // get number of chains
00298       int nchains = MCMCGetNChains();
00299 
00300       // get number of parameters
00301       int npar = GetNParameters();
00302       //cout << "number of parameters is: " << npar << endl;
00303 
00304       // loop over all chains and fill histogram
00305       for (int i = 0; i < nchains; ++i) {
00306          // get the current values of the parameters. These are
00307          // stored in fMCMCx.
00308 
00309          TIterator* setiter = fParamsPOI->createIterator();
00310          int j = 0;
00311 
00312          while(setiter->Next()){
00313 
00314             //check parameter names
00315             BCParameter * tempBCparam = GetParameter(j);
00316             //cout << "Chain " << i << "param: " << tempBCparam->GetName() << endl;
00317 
00318             const char * paramnamepointer = (tempBCparam->GetName()).c_str();
00319 
00320             double xij = fMCMCx.at(i * npar + j);
00321             AddToCurrentChainElement(xij, i, j);
00322             RooRealVar* parampointer = (RooRealVar*) &(_parametersForMarkovChainCurrent[paramnamepointer]);
00323             parampointer->setVal(xij);
00324             j++;
00325          }
00326          // will only work if _parametersForMarkovChain had correct info!
00327 
00328 //test stuff begin
00329 //var1->setVal( RooRandom::randomGenerator()->Uniform(100.) );
00330 //
00331 //_roostatsMarkovChain->Add(_parametersForMarkovChain_test, 0.001, 1.0);
00332 //
00333 //test stuff end
00334 
00335          if( !(EqualsLastChainElement(i)) ) {
00336             double weight = GetWeightForChain(i);
00337             _roostatsMarkovChain->Add(_parametersForMarkovChainPrevious, -1.* MCMCGetLogProbx(j), weight); 
00338             _parametersForMarkovChainPrevious = _parametersForMarkovChainCurrent;
00339          }
00340       }
00341    }
00342 }
00343 
00344 // ---------------------------------------------------------
00345 void BCRooInterface::AddToCurrentChainElement(double xij, int chainNum, int poiNum)
00346 {
00347    fCurrentStep[chainNum][poiNum] = xij;
00348 }
00349 
00350 // ---------------------------------------------------------
00351 bool BCRooInterface::EqualsLastChainElement(int chainNum)
00352 {
00353    bool equals = true;
00354    vector<double>::iterator itPrevious = fPreviousStep[chainNum].begin();
00355 
00356    if(fFirstComparison == true) {
00357       fFirstComparison = false;
00358       _parametersForMarkovChainPrevious = _parametersForMarkovChainCurrent;
00359       return true;
00360    }
00361 
00362 
00363    for (vector<double>::iterator itCurrent = fCurrentStep[chainNum].begin(); itCurrent != fCurrentStep[chainNum].end(); ++itCurrent) {
00364       if(*itCurrent != *itPrevious) {
00365          equals = false;
00366          fPreviousStep[chainNum] = fCurrentStep[chainNum];
00367          break;
00368       }
00369       ++itPrevious;
00370    }
00371 
00372    if(equals == true){
00373       fVecWeights[chainNum] += 1.0;
00374    }
00375 
00376    return equals;
00377 
00378 }
00379 
00380 // ---------------------------------------------------------
00381 double BCRooInterface::GetWeightForChain(int chainNum)
00382 {
00383    double retval = fVecWeights[chainNum];
00384    fVecWeights[chainNum]= 1.0 ;
00385    return retval;
00386 }
00387 

Generated by  doxygen 1.7.1