00001 #include "RooGlobalFunc.h"
00002 #include "RooMsgService.h"
00003 #include "RooProdPdf.h"
00004 #include "RooRealVar.h"
00005 #include "RooWorkspace.h"
00006 #include "TFile.h"
00007
00008 #include "BCRooInterface.h"
00009 #include <BAT/BCMath.h>
00010 #include <iostream>
00011
00012
00013
00014 void BCRooInterface::Initialize( const char* rootFile,
00015 const char* wsName,
00016 const char* dataName,
00017 const char* modelName,
00018 const char* priorName,
00019 const char* priorNuisanceName,
00020 const char* paramsName,
00021 const char* listPOIName )
00022 {
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 std::cout << "Opening " << rootFile << std::endl;
00038 TFile* file = new TFile(rootFile);
00039 std::cout << "content :\n";
00040 file->ls();
00041
00042 RooWorkspace* bat_ws = (RooWorkspace*) file->Get(wsName);
00043 bat_ws->Print("v");
00044
00045 fData = (RooAbsData*) bat_ws->data(dataName);
00046 fModel = (RooAbsPdf*) bat_ws->function(modelName);
00047
00048
00049 RooAbsPdf* priorPOI = (RooAbsPdf*) bat_ws->function(priorName);
00050 RooAbsPdf* priorNuisance = (RooAbsPdf*) bat_ws->pdf(priorNuisanceName);
00051 if (priorNuisance!=0 && priorPOI!=0) {
00052 fPrior = new RooProdPdf("fPrior","complete prior",*priorPOI,*priorNuisance);
00053 } else {
00054 if ( priorNuisance!=0 ) fPrior=priorNuisance;
00055 else if ( priorPOI!=0 ) fPrior = priorPOI;
00056 else std::cout << "No prior PDF: the program will crash\n";
00057 }
00058
00059 std::cout << "Imported parameters:\n";
00060 fParams = new RooArgList(*(bat_ws->set(listPOIName)));
00061 RooArgSet* paramsTmp = (RooArgSet*) bat_ws->set(paramsName);
00062 if (paramsTmp!=0) fParams->add(*paramsTmp);
00063 fParams->Print("v");
00064
00065
00066 fNll = new RooNLLVar("fNll","",*fModel,*fData,true);
00067
00068 file->Close();
00069
00070 DefineParameters();
00071 }
00072
00073
00074 BCRooInterface::BCRooInterface() : BCModel()
00075 {
00076
00077 }
00078
00079
00080 BCRooInterface::BCRooInterface(const char* name) : BCModel(name)
00081 {
00082
00083 }
00084
00085
00086 BCRooInterface::~BCRooInterface()
00087 {
00088
00089 }
00090
00091
00092 void BCRooInterface::DefineParameters()
00093 {
00094
00095 int default_nbins = 500;
00096
00097 int nParams = fParams->getSize();
00098 for (int iParam=0; iParam<nParams; iParam++) {
00099 RooRealVar* ipar = (RooRealVar*) fParams->at(iParam);
00100 this->AddParameter(ipar->GetName(),ipar->getMin(),ipar->getMax());
00101 this->SetNbins(ipar->GetName(),default_nbins);
00102 std::cout << "added parameter: " << ipar->GetName() << " defined in range [ " << ipar->getMin() << " - " << ipar->getMax() << " ]\n";
00103 }
00104 }
00105
00106
00107 double BCRooInterface::LogLikelihood(std::vector <double> parameters)
00108 {
00109
00110
00111 int nParams = fParams->getSize();
00112 for (int iParam=0; iParam<nParams; iParam++) {
00113 RooRealVar* ipar = (RooRealVar*) fParams->at(iParam);
00114 ipar->setVal(parameters.at(iParam));
00115 }
00116
00117
00118 double logprob = -fNll->getVal();
00119 return logprob;
00120 }
00121
00122
00123 double BCRooInterface::LogAPrioriProbability(std::vector <double> parameters)
00124 {
00125
00126
00127 int nParams = fParams->getSize();
00128 for (int iParam=0; iParam<nParams; iParam++) {
00129 RooRealVar* ipar = (RooRealVar*) fParams->at(iParam);
00130 ipar->setVal(parameters.at(iParam));
00131 }
00132
00133
00134 RooArgSet* tmpArgSet = new RooArgSet(*fParams);
00135 double prob = fPrior->getVal(tmpArgSet);
00136 delete tmpArgSet;
00137 if (prob<1e-300) prob = 1e-300;
00138 return log(prob);
00139 }
00140