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
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
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
00056
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
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
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
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
00135
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 {
00148 _default_nbins = 500;
00149 _fillChain = false;
00150 _addeddummyprior = false;
00151 }
00152
00153
00154 BCRooInterface::BCRooInterface(const char* name, bool fillChain) : BCModel(name)
00155 {
00156 _default_nbins = 500;
00157 _fillChain = fillChain;
00158 }
00159
00160
00161 BCRooInterface::~BCRooInterface()
00162 {
00163 if(_fillChain){
00164 delete _roostatsMarkovChain;
00165 }
00166
00167
00168
00169 }
00170
00171
00172 void BCRooInterface::DefineParameters()
00173 {
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 {
00193
00194
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
00202 double logprob = -fNll->getVal();
00203 return logprob;
00204 }
00205
00206
00207 double BCRooInterface::LogAPrioriProbability(const std::vector<double> & parameters)
00208 {
00209
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
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
00246
00247
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
00254
00255 _roostatsMarkovChain = new RooStats::MarkovChain();
00256
00257
00258
00259
00260 cout << "setting up parameters for RooStats markov chain" << endl;
00261 _parametersForMarkovChainPrevious.writeToStream(cout, false);
00262
00263
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
00281
00282
00283
00284
00285
00286
00287 }
00288
00289
00290
00291 void BCRooInterface::MCMCIterationInterface()
00292 {
00293
00294
00295 if(_fillChain){
00296
00297
00298 int nchains = MCMCGetNChains();
00299
00300
00301 int npar = GetNParameters();
00302
00303
00304
00305 for (int i = 0; i < nchains; ++i) {
00306
00307
00308
00309 TIterator* setiter = fParamsPOI->createIterator();
00310 int j = 0;
00311
00312 while(setiter->Next()){
00313
00314
00315 BCParameter * tempBCparam = GetParameter(j);
00316
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
00327
00328
00329
00330
00331
00332
00333
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