00001 #include "../interface/VerticalInterpHistPdf.h"
00002
00003 #include <cassert>
00004 #include <memory>
00005
00006 #include "RooFit.h"
00007 #include "Riostream.h"
00008
00009 #include "TIterator.h"
00010 #include "RooRealVar.h"
00011 #include "RooMsgService.h"
00012
00013
00014 #ifdef TRACE_CALLS
00015 #include "../interface/ProfilingTools.h"
00016 #define TRACEME() PerfCounter::add( __PRETTY_FUNCTION__ );
00017 #else
00018 #define TRACEME()
00019 #endif
00020
00021 ClassImp(VerticalInterpHistPdf)
00022
00023
00024
00025 VerticalInterpHistPdf::VerticalInterpHistPdf() :
00026 _cacheTotal(0),
00027 _cacheSingle(0),
00028 _cacheSingleGood(0)
00029 {
00030
00031 _funcIter = _funcList.createIterator() ;
00032 _coefIter = _coefList.createIterator() ;
00033 }
00034
00035
00036
00037 VerticalInterpHistPdf::VerticalInterpHistPdf(const char *name, const char *title, const RooRealVar &x, const RooArgList& inFuncList, const RooArgList& inCoefList, Double_t smoothRegion, Int_t smoothAlgo) :
00038 RooAbsPdf(name,title),
00039 _x("x","Independent variable",this,const_cast<RooRealVar&>(x)),
00040 _funcList("funcList","List of pdfs",this),
00041 _coefList("coefList","List of coefficients",this),
00042 _smoothRegion(smoothRegion),
00043 _smoothAlgo(smoothAlgo),
00044 _cacheTotal(0),
00045 _cacheSingle(0),
00046 _cacheSingleGood(0)
00047 {
00048
00049 if (inFuncList.getSize()!=2*inCoefList.getSize()+1) {
00050 coutE(InputArguments) << "VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName()
00051 << ") number of pdfs and coefficients inconsistent, must have Nfunc=1+2*Ncoef" << endl ;
00052 assert(0);
00053 }
00054
00055 TIterator* funcIter = inFuncList.createIterator() ;
00056 RooAbsArg* func;
00057 while((func = (RooAbsArg*)funcIter->Next())) {
00058 RooAbsPdf *pdf = dynamic_cast<RooAbsPdf*>(func);
00059 if (!pdf) {
00060 coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") function " << func->GetName() << " is not of type RooAbsPdf" << endl;
00061 assert(0);
00062 }
00063 RooArgSet *params = pdf->getParameters(RooArgSet(x));
00064 if (params->getSize() > 0) {
00065 coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") pdf " << func->GetName() << " has some parameters." << endl;
00066 assert(0);
00067 }
00068 delete params;
00069 _funcList.add(*func) ;
00070 }
00071 delete funcIter;
00072
00073 TIterator* coefIter = inCoefList.createIterator() ;
00074 RooAbsArg* coef;
00075 while((coef = (RooAbsArg*)coefIter->Next())) {
00076 if (!dynamic_cast<RooAbsReal*>(coef)) {
00077 coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal" << endl;
00078 assert(0);
00079 }
00080 _coefList.add(*coef) ;
00081 }
00082 delete coefIter;
00083
00084 _funcIter = _funcList.createIterator() ;
00085 _coefIter = _coefList.createIterator() ;
00086
00087 }
00088
00089
00090
00091
00092
00093 VerticalInterpHistPdf::VerticalInterpHistPdf(const VerticalInterpHistPdf& other, const char* name) :
00094 RooAbsPdf(other,name),
00095 _x("x",this,other._x),
00096 _funcList("funcList",this,other._funcList),
00097 _coefList("coefList",this,other._coefList),
00098 _smoothRegion(other._smoothRegion),
00099 _smoothAlgo(other._smoothAlgo),
00100 _cacheTotal(0),
00101 _cacheSingle(0),
00102 _cacheSingleGood(0)
00103 {
00104
00105
00106 _funcIter = _funcList.createIterator() ;
00107 _coefIter = _coefList.createIterator() ;
00108 }
00109
00110
00111
00112
00113 VerticalInterpHistPdf::~VerticalInterpHistPdf()
00114 {
00115
00116 delete _funcIter ;
00117 delete _coefIter ;
00118 if (_cacheTotal) {
00119 delete _cacheTotal;
00120 for (int i = 0; i < _funcList.getSize(); ++i) delete _cacheSingle[i];
00121 delete [] _cacheSingle;
00122 delete [] _cacheSingleGood;
00123 }
00124 }
00125
00126
00127 Double_t VerticalInterpHistPdf::evaluate() const
00128 {
00129 if (_cacheTotal == 0) setupCaches();
00130 #if 0
00131 printf("Evaluate called for x %f, cache status %d: ", _x.arg().getVal(), _sentry.good());
00132 int ndim = _coefList.getSize();
00133 for (int i = 0; i < ndim; ++i) {
00134 RooAbsReal *rar = dynamic_cast<RooAbsReal *>(_coefList.at(i));
00135 printf("%s = %+6.4f ", rar->GetName(), rar->getVal());
00136 }
00137 std::cout << std::endl;
00138 #endif
00139
00140 if (!_sentry.good()) syncTotal();
00141 int nbin = _cacheTotal->GetNbinsX(), ibin = _cacheTotal->FindBin(_x);
00142 if (ibin < 1) ibin = 1;
00143 else if (ibin > nbin) ibin = nbin;
00144 return _cacheTotal->GetBinContent(ibin);
00145 }
00146
00147
00148 void VerticalInterpHistPdf::syncComponent(int i) const {
00149 RooAbsPdf *pdfi = dynamic_cast<RooAbsPdf *>(_funcList.at(i));
00150 if (_cacheSingle[i] != 0) delete _cacheSingle[i];
00151 _cacheSingle[i] = pdfi->createHistogram("",dynamic_cast<const RooRealVar &>(_x.arg()));
00152 _cacheSingle[i]->SetDirectory(0);
00153 if (_cacheSingle[i]->Integral("width")) { _cacheSingle[i]->Scale(1.0/_cacheSingle[i]->Integral("width")); }
00154 if (i > 0) {
00155 for (int b = 1, nb = _cacheSingle[i]->GetNbinsX(); b <= nb; ++b) {
00156 double y = _cacheSingle[i]->GetBinContent(b);
00157 double y0 = _cacheSingle[0]->GetBinContent(b);
00158 if (_smoothAlgo < 0) {
00159 if (y > 0 && y0 > 0) {
00160 double logk = log(y/y0);
00161
00162
00163 _cacheSingle[i]->SetBinContent(b, logk);
00164 } else {
00165 _cacheSingle[i]->SetBinContent(b, 0);
00166 }
00167 } else {
00168 _cacheSingle[i]->SetBinContent(b, y - y0);
00169 }
00170 }
00171 }
00172 _cacheSingleGood[i] = true;
00173 }
00174
00175 void VerticalInterpHistPdf::syncTotal() const {
00176 int ndim = _coefList.getSize();
00177 for (int i = 0; i < 2*ndim+1; ++i) {
00178 if (!_cacheSingleGood[i]) syncComponent(i);
00179 }
00180 for (int b = 1, nb = _cacheTotal->GetNbinsX(); b <= nb; ++b) {
00181 double val = _cacheSingle[0]->GetBinContent(b);
00182 _coefIter->Reset();
00183 for (int i = 0; i < ndim; ++i) {
00184 double dhi = _cacheSingle[2*i+1]->GetBinContent(b);
00185 double dlo = _cacheSingle[2*i+2]->GetBinContent(b);
00186 double x = (dynamic_cast<RooAbsReal *>(_coefIter->Next()))->getVal();
00187 double alpha = x * 0.5 * ((dhi-dlo) + (dhi+dlo)*smoothStepFunc(x));
00188
00189
00190
00191
00192
00193
00194 if (_smoothAlgo < 0) {
00195 val *= exp(alpha);
00196 } else {
00197 val += alpha;
00198 }
00199 }
00200 if (val <= 0) val = 1e-9;
00201 _cacheTotal->SetBinContent(b, val);
00202 }
00203 double norm = _cacheTotal->Integral("width");
00204 if (norm > 0) _cacheTotal->Scale(1.0/norm);
00205 _sentry.reset();
00206 }
00207
00208 void VerticalInterpHistPdf::setupCaches() const {
00209 int ndim = _coefList.getSize();
00210 _cacheTotal = (dynamic_cast<const RooRealVar &>(_x.arg())).createHistogram("total");
00211 _cacheTotal->SetDirectory(0);
00212 _cacheSingle = new TH1*[2*ndim+1]; assert(_cacheSingle);
00213 _cacheSingleGood = new int[2*ndim+1]; assert(_cacheSingleGood);
00214 for (int i = 0; i < 2*ndim+1; ++i) {
00215 _cacheSingle[i] = 0;
00216 _cacheSingleGood[i] = 0;
00217 syncComponent(i);
00218 }
00219 if (_sentry.empty()) _sentry.addVars(_coefList);
00220 syncTotal();
00221 }
00222
00223
00224 ClassImp(FastVerticalInterpHistPdfBase)
00225 ClassImp(FastVerticalInterpHistPdf)
00226 ClassImp(FastVerticalInterpHistPdf2D)
00227
00228
00229
00230 FastVerticalInterpHistPdfBase::FastVerticalInterpHistPdfBase()
00231 {
00232
00233 _funcIter = _funcList.createIterator() ;
00234 _coefIter = _coefList.createIterator() ;
00235 }
00236
00237
00238
00239 FastVerticalInterpHistPdfBase::FastVerticalInterpHistPdfBase(const char *name, const char *title, const RooArgSet &obs, const RooArgList& inFuncList, const RooArgList& inCoefList, Double_t smoothRegion, Int_t smoothAlgo) :
00240 RooAbsPdf(name,title),
00241 _funcList("funcList","List of pdfs",this),
00242 _coefList("coefList","List of coefficients",this),
00243 _smoothRegion(smoothRegion),
00244 _smoothAlgo(smoothAlgo),
00245 _init(false),
00246 _morphs(), _morphParams()
00247 {
00248 TRACEME()
00249
00250 if (inFuncList.getSize()!=2*inCoefList.getSize()+1) {
00251 coutE(InputArguments) << "VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName()
00252 << ") number of pdfs and coefficients inconsistent, must have Nfunc=1+2*Ncoef" << endl ;
00253 assert(0);
00254 }
00255
00256 TIterator* funcIter = inFuncList.createIterator() ;
00257 RooAbsArg* func;
00258 while((func = (RooAbsArg*)funcIter->Next())) {
00259 RooAbsPdf *pdf = dynamic_cast<RooAbsPdf*>(func);
00260 if (!pdf) {
00261 coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") function " << func->GetName() << " is not of type RooAbsPdf" << endl;
00262 assert(0);
00263 }
00264 RooArgSet *params = pdf->getParameters(obs);
00265 if (params->getSize() > 0) {
00266 coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") pdf " << func->GetName() << " has some parameters." << endl;
00267 assert(0);
00268 }
00269 delete params;
00270 _funcList.add(*func) ;
00271 }
00272 delete funcIter;
00273
00274 TIterator* coefIter = inCoefList.createIterator() ;
00275 RooAbsArg* coef;
00276 while((coef = (RooAbsArg*)coefIter->Next())) {
00277 if (!dynamic_cast<RooAbsReal*>(coef)) {
00278 coutE(InputArguments) << "ERROR: VerticalInterpHistPdf::VerticalInterpHistPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal" << endl;
00279 assert(0);
00280 }
00281 _coefList.add(*coef) ;
00282 }
00283 delete coefIter;
00284
00285 _funcIter = _funcList.createIterator() ;
00286 _coefIter = _coefList.createIterator() ;
00287
00288 }
00289
00290
00291 FastVerticalInterpHistPdfBase::FastVerticalInterpHistPdfBase(const FastVerticalInterpHistPdfBase& other, const char* name) :
00292 RooAbsPdf(other,name),
00293 _funcList("funcList",this,other._funcList),
00294 _coefList("coefList",this,other._coefList),
00295 _smoothRegion(other._smoothRegion),
00296 _smoothAlgo(other._smoothAlgo),
00297 _init(false),
00298 #ifdef FastVerticalInterpHistPdf_CopyConstructor
00299 _morphs(other._morphs), _morphParams(other._morphParams)
00300 #else
00301 _morphs(), _morphParams()
00302 #endif
00303 {
00304
00305 _funcIter = _funcList.createIterator() ;
00306 _coefIter = _coefList.createIterator() ;
00307 #ifdef FastVerticalInterpHistPdf_CopyConstructor
00308 _sentry.addVars(_coefList);
00309 _sentry.setValueDirty();
00310 #endif
00311 }
00312
00313
00314
00315 FastVerticalInterpHistPdfBase::~FastVerticalInterpHistPdfBase()
00316 {
00317
00318 delete _funcIter ;
00319 delete _coefIter ;
00320 }
00321
00322
00323
00324 Double_t FastVerticalInterpHistPdf::evaluate() const
00325 {
00326 TRACEME()
00327 if (_cache.size() == 0) setupCaches();
00328 if (!_sentry.good() || !_init) syncTotal();
00329
00330 return _cache.GetAt(_x);
00331 }
00332
00333
00334 Double_t FastVerticalInterpHistPdf2D::evaluate() const
00335 {
00336 TRACEME()
00337 if (_cache.size() == 0) setupCaches();
00338
00339 if (!_sentry.good() || !_init) syncTotal();
00340
00341 return _cache.GetAt(_x, _y);
00342 }
00343
00344
00345
00346 void FastVerticalInterpHistPdf::syncNominal() const {
00347 TRACEME()
00348 RooAbsPdf *pdf = dynamic_cast<RooAbsPdf *>(_funcList.at(0));
00349 const RooRealVar &x = dynamic_cast<const RooRealVar &>(_x.arg());
00350 std::auto_ptr<TH1> hist(pdf->createHistogram("",x));
00351 hist->SetDirectory(0);
00352 _cacheNominal = FastHisto(*hist);
00353 _cacheNominal.Normalize();
00354 if (_smoothAlgo < 0) {
00355 _cacheNominalLog = _cacheNominal;
00356 _cacheNominalLog.Log();
00357 }
00358 }
00359
00360 void FastVerticalInterpHistPdf2D::syncNominal() const {
00361 TRACEME()
00362 RooAbsPdf *pdf = dynamic_cast<RooAbsPdf *>(_funcList.at(0));
00363 const RooRealVar &x = dynamic_cast<const RooRealVar &>(_x.arg());
00364 const RooRealVar &y = dynamic_cast<const RooRealVar &>(_y.arg());
00365 const RooCmdArg &cond = _conditional ? RooFit::ConditionalObservables(RooArgSet(x)) : RooCmdArg::none();
00366 std::auto_ptr<TH1> hist(pdf->createHistogram("", x, RooFit::YVar(y), cond));
00367 hist->SetDirectory(0);
00368 _cacheNominal = FastHisto2D(dynamic_cast<TH2F&>(*hist), _conditional);
00369 if (_conditional) _cacheNominal.NormalizeXSlices();
00370 else _cacheNominal.Normalize();
00371
00372 if (_smoothAlgo < 0) {
00373 _cacheNominalLog = _cacheNominal;
00374 _cacheNominalLog.Log();
00375 }
00376 }
00377
00378
00379
00380 void FastVerticalInterpHistPdfBase::syncMorph(Morph &out, const FastTemplate &nominal, FastTemplate &lo, FastTemplate &hi) const {
00381 if (_smoothAlgo < 0) {
00382 hi.LogRatio(nominal); lo.LogRatio(nominal);
00383
00384 } else {
00385 hi.Subtract(nominal); lo.Subtract(nominal);
00386
00387 }
00388 FastTemplate::SumDiff(hi, lo, out.sum, out.diff);
00389
00390 }
00391
00392 void FastVerticalInterpHistPdf::syncComponents(int dim) const {
00393 TRACEME()
00394 RooAbsPdf *pdfHi = dynamic_cast<RooAbsPdf *>(_funcList.at(2*dim+1));
00395 RooAbsPdf *pdfLo = dynamic_cast<RooAbsPdf *>(_funcList.at(2*dim+2));
00396 const RooRealVar &x = dynamic_cast<const RooRealVar &>(_x.arg());
00397 std::auto_ptr<TH1> histHi(pdfHi->createHistogram("",x)); histHi->SetDirectory(0);
00398 std::auto_ptr<TH1> histLo(pdfLo->createHistogram("",x)); histLo->SetDirectory(0);
00399
00400 FastHisto hi(*histHi), lo(*histLo);
00401
00402 hi.Normalize(); lo.Normalize();
00403
00404 syncMorph(_morphs[dim], _cacheNominal, lo, hi);
00405 }
00406 void FastVerticalInterpHistPdf2D::syncComponents(int dim) const {
00407 RooAbsPdf *pdfHi = dynamic_cast<RooAbsPdf *>(_funcList.at(2*dim+1));
00408 RooAbsPdf *pdfLo = dynamic_cast<RooAbsPdf *>(_funcList.at(2*dim+2));
00409 const RooRealVar &x = dynamic_cast<const RooRealVar &>(_x.arg());
00410 const RooRealVar &y = dynamic_cast<const RooRealVar &>(_y.arg());
00411 const RooCmdArg &cond = _conditional ? RooFit::ConditionalObservables(RooArgSet(x)) : RooCmdArg::none();
00412 std::auto_ptr<TH1> histHi(pdfHi->createHistogram("", x, RooFit::YVar(y), cond)); histHi->SetDirectory(0);
00413 std::auto_ptr<TH1> histLo(pdfLo->createHistogram("", x, RooFit::YVar(y), cond)); histLo->SetDirectory(0);
00414
00415 FastHisto2D hi(dynamic_cast<TH2&>(*histHi), _conditional), lo(dynamic_cast<TH2&>(*histLo), _conditional);
00416
00417 if (_conditional) {
00418 hi.NormalizeXSlices(); lo.NormalizeXSlices();
00419 } else {
00420 hi.Normalize(); lo.Normalize();
00421 }
00422
00423 syncMorph(_morphs[dim], _cacheNominal, lo, hi);
00424 }
00425
00426
00427 void FastVerticalInterpHistPdfBase::syncTotal(FastTemplate &cache, const FastTemplate &cacheNominal, const FastTemplate &cacheNominalLog) const {
00428 TRACEME()
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447 cache.CopyValues(_smoothAlgo < 0 ? cacheNominalLog : cacheNominal);
00448
00449
00450
00451 for (int i = 0, ndim = _coefList.getSize(); i < ndim; ++i) {
00452 double x = _morphParams[i]->getVal();
00453 double a = 0.5*x, b = smoothStepFunc(x);
00454 cache.Meld(_morphs[i].diff, _morphs[i].sum, a, b);
00455
00456 }
00457
00458
00459 if (_smoothAlgo < 0) {
00460 cache.Exp();
00461
00462 } else {
00463 cache.CropUnderflows();
00464 }
00465
00466
00467 _sentry.reset();
00468 _init = true;
00469 }
00470
00471 void FastVerticalInterpHistPdf::syncTotal() const {
00472 FastVerticalInterpHistPdfBase::syncTotal(_cache, _cacheNominal, _cacheNominalLog);
00473
00474
00475 _cache.Normalize();
00476
00477 }
00478
00479 void FastVerticalInterpHistPdf2D::syncTotal() const {
00480 FastVerticalInterpHistPdfBase::syncTotal(_cache, _cacheNominal, _cacheNominalLog);
00481
00482 if (_conditional) _cache.NormalizeXSlices();
00483 else _cache.Normalize();
00484
00485 }
00486
00487
00488
00489 void FastVerticalInterpHistPdf::setupCaches() const {
00490 TRACEME()
00491 int ndim = _coefList.getSize();
00492
00493 _morphs.resize(ndim);
00494 _morphParams.resize(ndim);
00495 syncNominal();
00496
00497 _coefIter->Reset();
00498 for (int i = 0; i < ndim; ++i) {
00499 _morphParams[i] = dynamic_cast<RooAbsReal *>(_coefIter->Next());
00500 _morphs[i].sum.Resize(_cacheNominal.size());
00501 _morphs[i].diff.Resize(_cacheNominal.size());
00502 syncComponents(i);
00503 }
00504 _cache = FastHisto(_cacheNominal);
00505
00506 if (_sentry.empty()) _sentry.addVars(_coefList);
00507 syncTotal();
00508 }
00509
00510 void FastVerticalInterpHistPdf2D::setupCaches() const {
00511 TRACEME()
00512 int ndim = _coefList.getSize();
00513
00514 _morphs.resize(ndim);
00515 _morphParams.resize(ndim);
00516 syncNominal();
00517
00518 _coefIter->Reset();
00519 for (int i = 0; i < ndim; ++i) {
00520 _morphParams[i] = dynamic_cast<RooAbsReal *>(_coefIter->Next());
00521 _morphs[i].sum.Resize(_cacheNominal.size());
00522 _morphs[i].diff.Resize(_cacheNominal.size());
00523 syncComponents(i);
00524 }
00525 _cache = FastHisto2D(_cacheNominal);
00526
00527 if (_sentry.empty()) _sentry.addVars(_coefList);
00528 syncTotal();
00529 }
00530
00531