00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "RooFit.h"
00028 #include "Riostream.h"
00029
00030 #include "TH1.h"
00031 #include "TH1.h"
00032 #include "TDirectory.h"
00033 #include "TMath.h"
00034 #include "RooMsgService.h"
00035 #include "RooDataHist.h"
00036 #include "RooDataHistSliceIter.h"
00037 #include "RooAbsLValue.h"
00038 #include "RooArgList.h"
00039 #include "RooRealVar.h"
00040 #include "RooMath.h"
00041 #include "RooBinning.h"
00042 #include "RooPlot.h"
00043 #include "RooHistError.h"
00044 #include "RooCategory.h"
00045 #include "RooCmdConfig.h"
00046 #include "RooTreeDataStore.h"
00047 #include "TTree.h"
00048 #include "RooTreeData.h"
00049
00050 using namespace std ;
00051
00052 ClassImp(RooDataHist)
00053 ;
00054
00055
00056
00057
00058 RooDataHist::RooDataHist() : _pbinvCacheMgr(0,10)
00059 {
00060
00061 _arrSize = 0 ;
00062 _wgt = 0 ;
00063 _errLo = 0 ;
00064 _errHi = 0 ;
00065 _sumw2 = 0 ;
00066 _binv = 0 ;
00067 _pbinv = 0 ;
00068 _curWeight = 0 ;
00069 _curIndex = -1 ;
00070 _realIter = _realVars.createIterator() ;
00071 _binValid = 0 ;
00072
00073 }
00074
00075
00076
00077
00078 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const char* binningName) :
00079 RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00080 {
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 _dstore = new RooTreeDataStore(name,title,_vars) ;
00099
00100 initialize(binningName) ;
00101
00102 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00103
00104 appendToDir(this,kTRUE) ;
00105 }
00106
00107
00108
00109
00110 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const RooAbsData& data, Double_t wgt) :
00111 RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00112 {
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 _dstore = new RooTreeDataStore(name,title,_vars) ;
00135 initialize() ;
00136 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00137
00138 add(data,(const RooFormulaVar*)0,wgt) ;
00139 appendToDir(this,kTRUE) ;
00140 }
00141
00142
00143
00144
00145 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, RooCategory& indexCat,
00146 map<string,TH1*> histMap, Double_t wgt) :
00147 RooAbsData(name,title,RooArgSet(vars,&indexCat)),
00148 _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00149 {
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 _dstore = new RooTreeDataStore(name,title,_vars) ;
00160
00161 importTH1Set(vars, indexCat, histMap, wgt, kTRUE) ;
00162
00163 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00164 }
00165
00166
00167
00168
00169 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, RooCategory& indexCat,
00170 map<string,RooDataHist*> dhistMap, Double_t wgt) :
00171 RooAbsData(name,title,RooArgSet(vars,&indexCat)),
00172 _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00173 {
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 _dstore = new RooTreeDataStore(name,title,_vars) ;
00184
00185 importDHistSet(vars, indexCat, dhistMap, wgt) ;
00186
00187 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00188 }
00189
00190
00191
00192
00193 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const TH1* hist, Double_t wgt) :
00194 RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00195 {
00196
00197
00198
00199
00200
00201
00202 _dstore = new RooTreeDataStore(name,title,_vars) ;
00203
00204
00205 if (vars.getSize() != hist->GetDimension()) {
00206 coutE(InputArguments) << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
00207 << "number of dimension variables" << endl ;
00208 assert(0) ;
00209 }
00210
00211 importTH1(vars,*const_cast<TH1*>(hist),wgt, kTRUE) ;
00212
00213 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00214 }
00215
00216
00217
00218
00219 RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3,
00220 const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,const RooCmdArg& arg7,const RooCmdArg& arg8) :
00221 RooAbsData(name,title,RooArgSet(vars,(RooAbsArg*)RooCmdConfig::decodeObjOnTheFly("RooDataHist::RooDataHist", "IndexCat",0,0,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8))),
00222 _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00223 {
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 _dstore = new RooTreeDataStore(name,title,_vars) ;
00259
00260
00261 RooCmdConfig pc(Form("RooDataHist::ctor(%s)",GetName())) ;
00262 pc.defineObject("impHist","ImportHisto",0) ;
00263 pc.defineInt("impDens","ImportHisto",0) ;
00264 pc.defineObject("indexCat","IndexCat",0) ;
00265 pc.defineObject("impSliceHist","ImportHistoSlice",0,0,kTRUE) ;
00266 pc.defineString("impSliceState","ImportHistoSlice",0,"",kTRUE) ;
00267 pc.defineObject("impSliceDHist","ImportDataHistSlice",0,0,kTRUE) ;
00268 pc.defineString("impSliceDState","ImportDataHistSlice",0,"",kTRUE) ;
00269 pc.defineDouble("weight","Weight",0,1) ;
00270 pc.defineObject("dummy1","ImportDataHistSliceMany",0) ;
00271 pc.defineObject("dummy2","ImportHistoSliceMany",0) ;
00272 pc.defineMutex("ImportHisto","ImportHistoSlice","ImportDataHistSlice") ;
00273 pc.defineDependency("ImportHistoSlice","IndexCat") ;
00274 pc.defineDependency("ImportDataHistSlice","IndexCat") ;
00275
00276 RooLinkedList l ;
00277 l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
00278 l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
00279 l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
00280 l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
00281
00282
00283 pc.process(l) ;
00284 if (!pc.ok(kTRUE)) {
00285 assert(0) ;
00286 return ;
00287 }
00288
00289 TH1* impHist = static_cast<TH1*>(pc.getObject("impHist")) ;
00290 Bool_t impDens = pc.getInt("impDens") ;
00291 Double_t initWgt = pc.getDouble("weight") ;
00292 const char* impSliceNames = pc.getString("impSliceState","",kTRUE) ;
00293 const RooLinkedList& impSliceHistos = pc.getObjectList("impSliceHist") ;
00294 RooCategory* indexCat = static_cast<RooCategory*>(pc.getObject("indexCat")) ;
00295 const char* impSliceDNames = pc.getString("impSliceDState","",kTRUE) ;
00296 const RooLinkedList& impSliceDHistos = pc.getObjectList("impSliceDHist") ;
00297
00298
00299 if (impHist) {
00300
00301
00302 importTH1(vars,*impHist,initWgt, impDens) ;
00303
00304 } else if (indexCat) {
00305
00306
00307 if (impSliceHistos.GetSize()>0) {
00308
00309
00310 map<string,TH1*> hmap ;
00311 char tmp[1024] ;
00312 strlcpy(tmp,impSliceNames,1024) ;
00313 char* token = strtok(tmp,",") ;
00314 TIterator* hiter = impSliceHistos.MakeIterator() ;
00315 while(token) {
00316 hmap[token] = (TH1*) hiter->Next() ;
00317 token = strtok(0,",") ;
00318 }
00319 importTH1Set(vars,*indexCat,hmap,initWgt,kTRUE) ;
00320 } else {
00321
00322
00323 map<string,RooDataHist*> dmap ;
00324 char tmp[1024] ;
00325 strlcpy(tmp,impSliceDNames,1024) ;
00326 char* token = strtok(tmp,",") ;
00327 TIterator* hiter = impSliceDHistos.MakeIterator() ;
00328 while(token) {
00329 dmap[token] = (RooDataHist*) hiter->Next() ;
00330 token = strtok(0,",") ;
00331 }
00332 importDHistSet(vars,*indexCat,dmap,initWgt) ;
00333
00334 }
00335
00336
00337 } else {
00338
00339
00340 initialize() ;
00341 appendToDir(this,kTRUE) ;
00342
00343 }
00344
00345 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00346
00347 }
00348
00349
00350
00351
00352
00353 void RooDataHist::importTH1(const RooArgList& vars, TH1& histo, Double_t wgt, Bool_t doDensityCorrection)
00354 {
00355
00356
00357
00358 Int_t offset[3] ;
00359 adjustBinning(vars,histo,offset) ;
00360
00361
00362 initialize() ;
00363 appendToDir(this,kTRUE) ;
00364
00365
00366 RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
00367 RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
00368 RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
00369
00370
00371 Int_t xmin(0),ymin(0),zmin(0) ;
00372 RooArgSet vset(*xvar) ;
00373 Double_t volume = xvar->getMax()-xvar->getMin() ;
00374 xmin = offset[0] ;
00375 if (yvar) {
00376 vset.add(*yvar) ;
00377 ymin = offset[1] ;
00378 volume *= (yvar->getMax()-yvar->getMin()) ;
00379 }
00380 if (zvar) {
00381 vset.add(*zvar) ;
00382 zmin = offset[2] ;
00383 volume *= (zvar->getMax()-zvar->getMin()) ;
00384 }
00385 Double_t avgBV = volume / numEntries() ;
00386
00387 Int_t ix(0),iy(0),iz(0) ;
00388 for (ix=0 ; ix < xvar->getBins() ; ix++) {
00389 xvar->setBin(ix) ;
00390 if (yvar) {
00391 for (iy=0 ; iy < yvar->getBins() ; iy++) {
00392 yvar->setBin(iy) ;
00393 if (zvar) {
00394 for (iz=0 ; iz < zvar->getBins() ; iz++) {
00395 zvar->setBin(iz) ;
00396 Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
00397 add(vset,bv*histo.GetBinContent(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,2)) ;
00398 }
00399 } else {
00400 Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
00401 add(vset,bv*histo.GetBinContent(ix+1+xmin,iy+1+ymin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin,iy+1+ymin)*wgt,2)) ;
00402 }
00403 }
00404 } else {
00405 Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
00406 add(vset,bv*histo.GetBinContent(ix+1+xmin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin)*wgt,2)) ;
00407 }
00408 }
00409
00410 }
00411
00412
00413
00414
00415
00416
00417 void RooDataHist::importTH1Set(const RooArgList& vars, RooCategory& indexCat, map<string,TH1*> hmap, Double_t wgt, Bool_t doDensityCorrection)
00418 {
00419
00420
00421
00422
00423 RooCategory* icat = (RooCategory*) _vars.find(indexCat.GetName()) ;
00424
00425 TH1* histo(0) ;
00426 Bool_t init(kFALSE) ;
00427 for (map<string,TH1*>::iterator hiter = hmap.begin() ; hiter!=hmap.end() ; ++hiter) {
00428
00429 if (!histo) {
00430 histo = hiter->second ;
00431 }
00432
00433 if (!indexCat.lookupType(hiter->first.c_str())) {
00434 indexCat.defineType(hiter->first.c_str()) ;
00435 coutI(InputArguments) << "RooDataHist::importTH1Set(" << GetName() << ") defining state \"" << hiter->first << "\" in index category " << indexCat.GetName() << endl ;
00436 }
00437 if (!icat->lookupType(hiter->first.c_str())) {
00438 icat->defineType(hiter->first.c_str()) ;
00439 }
00440 }
00441
00442
00443 if (histo && (vars.getSize() != histo->GetDimension())) {
00444 coutE(InputArguments) << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
00445 << "number of continuous variables" << endl ;
00446 assert(0) ;
00447 }
00448
00449
00450 Int_t offset[3] ;
00451 adjustBinning(vars,*histo,offset) ;
00452
00453
00454 if (!init) {
00455 initialize() ;
00456 appendToDir(this,kTRUE) ;
00457 init = kTRUE ;
00458 }
00459
00460
00461 RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
00462 RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
00463 RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
00464
00465
00466 Int_t xmin(0),ymin(0),zmin(0) ;
00467 RooArgSet vset(*xvar) ;
00468 Double_t volume = xvar->getMax()-xvar->getMin() ;
00469 xmin = offset[0] ;
00470 if (yvar) {
00471 vset.add(*yvar) ;
00472 ymin = offset[1] ;
00473 volume *= (yvar->getMax()-yvar->getMin()) ;
00474 }
00475 if (zvar) {
00476 vset.add(*zvar) ;
00477 zmin = offset[2] ;
00478 volume *= (zvar->getMax()-zvar->getMin()) ;
00479 }
00480 Double_t avgBV = volume / numEntries() ;
00481
00482 Int_t ic(0),ix(0),iy(0),iz(0) ;
00483 for (ic=0 ; ic < icat->numBins(0) ; ic++) {
00484 icat->setBin(ic) ;
00485 histo = hmap[icat->getLabel()] ;
00486 for (ix=0 ; ix < xvar->getBins() ; ix++) {
00487 xvar->setBin(ix) ;
00488 if (yvar) {
00489 for (iy=0 ; iy < yvar->getBins() ; iy++) {
00490 yvar->setBin(iy) ;
00491 if (zvar) {
00492 for (iz=0 ; iz < zvar->getBins() ; iz++) {
00493 zvar->setBin(iz) ;
00494 Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
00495 add(vset,bv*histo->GetBinContent(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,2)) ;
00496 }
00497 } else {
00498 Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
00499 add(vset,bv*histo->GetBinContent(ix+1+xmin,iy+1+ymin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin,iy+1+ymin)*wgt,2)) ;
00500 }
00501 }
00502 } else {
00503 Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
00504 add(vset,bv*histo->GetBinContent(ix+1+xmin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin)*wgt,2)) ;
00505 }
00506 }
00507 }
00508
00509 }
00510
00511
00512
00513
00514 void RooDataHist::importDHistSet(const RooArgList& , RooCategory& indexCat, std::map<std::string,RooDataHist*> dmap, Double_t initWgt)
00515 {
00516
00517
00518
00519
00520 RooCategory* icat = (RooCategory*) _vars.find(indexCat.GetName()) ;
00521
00522 for (map<string,RooDataHist*>::iterator diter = dmap.begin() ; diter!=dmap.end() ; ++diter) {
00523
00524
00525 if (!indexCat.lookupType(diter->first.c_str())) {
00526 indexCat.defineType(diter->first.c_str()) ;
00527 coutI(InputArguments) << "RooDataHist::importDHistSet(" << GetName() << ") defining state \"" << diter->first << "\" in index category " << indexCat.GetName() << endl ;
00528 }
00529 if (!icat->lookupType(diter->first.c_str())) {
00530 icat->defineType(diter->first.c_str()) ;
00531 }
00532 }
00533
00534 initialize() ;
00535 appendToDir(this,kTRUE) ;
00536
00537
00538 for (map<string,RooDataHist*>::iterator diter = dmap.begin() ; diter!=dmap.end() ; ++diter) {
00539
00540 RooDataHist* dhist = diter->second ;
00541
00542 icat->setLabel(diter->first.c_str()) ;
00543
00544
00545 for (Int_t i=0 ; i<dhist->numEntries() ; i++) {
00546 _vars = *dhist->get(i) ;
00547 add(_vars,dhist->weight()*initWgt, pow(dhist->weightError(SumW2),2) ) ;
00548 }
00549
00550 }
00551 }
00552
00553
00554 void RooDataHist::adjustBinning(const RooArgList& vars, TH1& href, Int_t* offset)
00555 {
00556
00557
00558
00559
00560
00561 RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
00562 if (!dynamic_cast<RooRealVar*>(xvar)) {
00563 coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << xvar->GetName() << " must be real" << endl ;
00564 assert(0) ;
00565 }
00566
00567 Double_t xlo = ((RooRealVar*)vars.at(0))->getMin() ;
00568 Double_t xhi = ((RooRealVar*)vars.at(0))->getMax() ;
00569 Int_t xmin(0) ;
00570 if (href.GetXaxis()->GetXbins()->GetArray()) {
00571
00572 RooBinning xbins(href.GetNbinsX(),href.GetXaxis()->GetXbins()->GetArray()) ;
00573
00574
00575 Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+1e-6)) ;
00576 Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-1e-6)) ;
00577 xbins.setRange(xloAdj,xhiAdj) ;
00578 ((RooRealVar*)vars.at(0))->setBinning(xbins) ;
00579 if (fabs(xloAdj-xlo)>1e-6||fabs(xhiAdj-xhi)) {
00580 coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: ["
00581 << xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
00582 }
00583
00584 xvar->setBinning(xbins) ;
00585 xmin = xbins.rawBinNumber(xloAdj+1e-6) ;
00586 if (offset) {
00587 offset[0] = xmin ;
00588 }
00589
00590 } else {
00591 RooBinning xbins(href.GetXaxis()->GetXmin(),href.GetXaxis()->GetXmax()) ;
00592 xbins.addUniform(href.GetNbinsX(),href.GetXaxis()->GetXmin(),href.GetXaxis()->GetXmax()) ;
00593
00594
00595 Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+1e-6)) ;
00596 Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-1e-6)) ;
00597 xbins.setRange(xloAdj,xhiAdj) ;
00598 ((RooRealVar*)vars.at(0))->setRange(xloAdj,xhiAdj) ;
00599 if (fabs(xloAdj-xlo)>1e-6||fabs(xhiAdj-xhi)) {
00600 coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: ["
00601 << xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
00602 }
00603
00604
00605 RooUniformBinning xbins2(xloAdj,xhiAdj,xbins.numBins()) ;
00606 xvar->setBinning(xbins2) ;
00607 xmin = xbins.rawBinNumber(xloAdj+1e-6) ;
00608 if (offset) {
00609 offset[0] = xmin ;
00610 }
00611 }
00612
00613
00614
00615
00616 RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
00617 Int_t ymin(0) ;
00618 if (yvar) {
00619 Double_t ylo = ((RooRealVar*)vars.at(1))->getMin() ;
00620 Double_t yhi = ((RooRealVar*)vars.at(1))->getMax() ;
00621
00622 if (!dynamic_cast<RooRealVar*>(yvar)) {
00623 coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << yvar->GetName() << " must be real" << endl ;
00624 assert(0) ;
00625 }
00626
00627 if (href.GetYaxis()->GetXbins()->GetArray()) {
00628
00629 RooBinning ybins(href.GetNbinsY(),href.GetYaxis()->GetXbins()->GetArray()) ;
00630
00631
00632 Double_t yloAdj = ybins.binLow(ybins.binNumber(ylo+1e-6)) ;
00633 Double_t yhiAdj = ybins.binHigh(ybins.binNumber(yhi-1e-6)) ;
00634 ybins.setRange(yloAdj,yhiAdj) ;
00635 ((RooRealVar*)vars.at(1))->setBinning(ybins) ;
00636 if (fabs(yloAdj-ylo)>1e-6||fabs(yhiAdj-yhi)) {
00637 coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << yvar->GetName() << " expanded to nearest bin boundaries: ["
00638 << ylo << "," << yhi << "] --> [" << yloAdj << "," << yhiAdj << "]" << endl ;
00639 }
00640
00641 yvar->setBinning(ybins) ;
00642 ymin = ybins.rawBinNumber(yloAdj+1e-6) ;
00643 if (offset) {
00644 offset[1] = ymin ;
00645 }
00646
00647 } else {
00648
00649 RooBinning ybins(href.GetYaxis()->GetXmin(),href.GetYaxis()->GetXmax()) ;
00650 ybins.addUniform(href.GetNbinsY(),href.GetYaxis()->GetXmin(),href.GetYaxis()->GetXmax()) ;
00651
00652
00653 Double_t yloAdj = ybins.binLow(ybins.binNumber(ylo+1e-6)) ;
00654 Double_t yhiAdj = ybins.binHigh(ybins.binNumber(yhi-1e-6)) ;
00655 ybins.setRange(yloAdj,yhiAdj) ;
00656 ((RooRealVar*)vars.at(1))->setRange(yloAdj,yhiAdj) ;
00657 if (fabs(yloAdj-ylo)>1e-6||fabs(yhiAdj-yhi)) {
00658 coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << yvar->GetName() << " expanded to nearest bin boundaries: ["
00659 << ylo << "," << yhi << "] --> [" << yloAdj << "," << yhiAdj << "]" << endl ;
00660 }
00661
00662 RooUniformBinning ybins2(yloAdj,yhiAdj,ybins.numBins()) ;
00663 yvar->setBinning(ybins2) ;
00664 ymin = ybins.rawBinNumber(yloAdj+1e-6) ;
00665 if (offset) {
00666 offset[1] = ymin ;
00667 }
00668
00669 }
00670 }
00671
00672
00673 RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
00674 Int_t zmin(0) ;
00675 if (zvar) {
00676 Double_t zlo = ((RooRealVar*)vars.at(2))->getMin() ;
00677 Double_t zhi = ((RooRealVar*)vars.at(2))->getMax() ;
00678
00679 if (!dynamic_cast<RooRealVar*>(zvar)) {
00680 coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << zvar->GetName() << " must be real" << endl ;
00681 assert(0) ;
00682 }
00683
00684 if (href.GetZaxis()->GetXbins()->GetArray()) {
00685
00686 RooBinning zbins(href.GetNbinsZ(),href.GetZaxis()->GetXbins()->GetArray()) ;
00687
00688
00689 Double_t zloAdj = zbins.binLow(zbins.binNumber(zlo+1e-6)) ;
00690 Double_t zhiAdj = zbins.binHigh(zbins.binNumber(zhi-1e-6)) ;
00691 zbins.setRange(zloAdj,zhiAdj) ;
00692 ((RooRealVar*)vars.at(2))->setBinning(zbins) ;
00693 if (fabs(zloAdj-zlo)>1e-6||fabs(zhiAdj-zhi)) {
00694 coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << zvar->GetName() << " expanded to nearest bin boundaries: ["
00695 << zlo << "," << zhi << "] --> [" << zloAdj << "," << zhiAdj << "]" << endl ;
00696 }
00697
00698 zvar->setBinning(zbins) ;
00699 zmin = zbins.rawBinNumber(zloAdj+1e-6) ;
00700 if (offset) {
00701 offset[2] = zmin ;
00702 }
00703
00704 } else {
00705
00706 RooBinning zbins(href.GetZaxis()->GetXmin(),href.GetZaxis()->GetXmax()) ;
00707 zbins.addUniform(href.GetNbinsZ(),href.GetZaxis()->GetXmin(),href.GetZaxis()->GetXmax()) ;
00708
00709
00710 Double_t zloAdj = zbins.binLow(zbins.binNumber(zlo+1e-6)) ;
00711 Double_t zhiAdj = zbins.binHigh(zbins.binNumber(zhi-1e-6)) ;
00712 zbins.setRange(zloAdj,zhiAdj) ;
00713 ((RooRealVar*)vars.at(2))->setRange(zloAdj,zhiAdj) ;
00714 if (fabs(zloAdj-zlo)>1e-6||fabs(zhiAdj-zhi)) {
00715 coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << zvar->GetName() << " expanded to nearest bin boundaries: ["
00716 << zlo << "," << zhi << "] --> [" << zloAdj << "," << zhiAdj << "]" << endl ;
00717 }
00718
00719 RooUniformBinning zbins2(zloAdj,zhiAdj,zbins.numBins()) ;
00720 zvar->setBinning(zbins2) ;
00721 zmin = zbins.rawBinNumber(zloAdj+1e-6) ;
00722 if (offset) {
00723 offset[2] = zmin ;
00724 }
00725 }
00726 }
00727
00728 }
00729
00730
00731
00732
00733
00734
00735 void RooDataHist::initialize(const char* binningName, Bool_t fillTree)
00736 {
00737
00738
00739
00740
00741
00742
00743 RooAbsArg* real ;
00744 _iterator->Reset() ;
00745 while((real=(RooAbsArg*)_iterator->Next())) {
00746 if (dynamic_cast<RooAbsReal*>(real)) _realVars.add(*real) ;
00747 }
00748 _realIter = _realVars.createIterator() ;
00749
00750
00751 _iterator->Reset() ;
00752 RooAbsArg* rvarg ;
00753 while((rvarg=(RooAbsArg*)_iterator->Next())) {
00754 if (binningName) {
00755 RooRealVar* rrv = dynamic_cast<RooRealVar*>(rvarg) ;
00756 if (rrv) {
00757 rrv->setBinning(rrv->getBinning(binningName)) ;
00758 }
00759 }
00760
00761 _lvvars.push_back(dynamic_cast<RooAbsLValue*>(rvarg)) ;
00762
00763 const RooAbsBinning* binning = dynamic_cast<RooAbsLValue*>(rvarg)->getBinningPtr(0) ;
00764 _lvbins.push_back(binning ? binning->clone() : 0) ;
00765 }
00766
00767
00768
00769 _idxMult.resize(_vars.getSize()) ;
00770
00771 _arrSize = 1 ;
00772 _iterator->Reset() ;
00773 RooAbsLValue* arg ;
00774 Int_t n(0), i ;
00775 while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
00776
00777
00778 for (i=0 ; i<n ; i++) {
00779 _idxMult[i] *= arg->numBins() ;
00780 }
00781 _idxMult[n++] = 1 ;
00782
00783
00784 _arrSize *= arg->numBins() ;
00785 }
00786
00787
00788 if (!_wgt) {
00789 _wgt = new Double_t[_arrSize] ;
00790 _errLo = new Double_t[_arrSize] ;
00791 _errHi = new Double_t[_arrSize] ;
00792 _sumw2 = new Double_t[_arrSize] ;
00793 _binv = new Double_t[_arrSize] ;
00794
00795
00796
00797 if (fillTree==kFALSE) {
00798 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00799 }
00800
00801 for (i=0 ; i<_arrSize ; i++) {
00802 _wgt[i] = 0 ;
00803 _errLo[i] = -1 ;
00804 _errHi[i] = -1 ;
00805 _sumw2[i] = 0 ;
00806 }
00807 }
00808
00809 if (!fillTree) return ;
00810
00811
00812
00813
00814 Int_t ibin ;
00815 for (ibin=0 ; ibin<_arrSize ; ibin++) {
00816 _iterator->Reset() ;
00817 RooAbsLValue* arg2 ;
00818 Int_t j(0), idx(0), tmp(ibin) ;
00819 Double_t theBinVolume(1) ;
00820 while((arg2=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
00821 idx = tmp / _idxMult[j] ;
00822 tmp -= idx*_idxMult[j++] ;
00823 RooAbsLValue* arglv = dynamic_cast<RooAbsLValue*>(arg2) ;
00824 arglv->setBin(idx) ;
00825 theBinVolume *= arglv->getBinWidth(idx) ;
00826 }
00827 _binv[ibin] = theBinVolume ;
00828 fill() ;
00829 }
00830
00831
00832 }
00833
00834
00835
00836
00837 RooDataHist::RooDataHist(const RooDataHist& other, const char* newname) :
00838 RooAbsData(other,newname), RooDirItem(), _idxMult(other._idxMult), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(other._pbinvCacheMgr,0)
00839 {
00840
00841
00842 Int_t i ;
00843
00844
00845 _arrSize = other._arrSize ;
00846 _wgt = new Double_t[_arrSize] ;
00847 _errLo = new Double_t[_arrSize] ;
00848 _errHi = new Double_t[_arrSize] ;
00849 _binv = new Double_t[_arrSize] ;
00850 _sumw2 = new Double_t[_arrSize] ;
00851 for (i=0 ; i<_arrSize ; i++) {
00852 _wgt[i] = other._wgt[i] ;
00853 _errLo[i] = other._errLo[i] ;
00854 _errHi[i] = other._errHi[i] ;
00855 _sumw2[i] = other._sumw2[i] ;
00856 _binv[i] = other._binv[i] ;
00857 }
00858
00859
00860 RooAbsArg* arg ;
00861 _iterator->Reset() ;
00862 while((arg=(RooAbsArg*)_iterator->Next())) {
00863 if (dynamic_cast<RooAbsReal*>(arg)) _realVars.add(*arg) ;
00864 }
00865 _realIter = _realVars.createIterator() ;
00866
00867
00868 _iterator->Reset() ;
00869 RooAbsArg* rvarg ;
00870 while((rvarg=(RooAbsArg*)_iterator->Next())) {
00871
00872 _lvvars.push_back(dynamic_cast<RooAbsLValue*>(rvarg)) ;
00873
00874 const RooAbsBinning* binning = dynamic_cast<RooAbsLValue*>(rvarg)->getBinningPtr(0) ;
00875 _lvbins.push_back(binning ? binning->clone() : 0) ;
00876 }
00877
00878 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00879
00880 appendToDir(this,kTRUE) ;
00881 }
00882
00883
00884
00885
00886 RooDataHist::RooDataHist(const char* name, const char* title, RooDataHist* h, const RooArgSet& varSubset,
00887 const RooFormulaVar* cutVar, const char* cutRange, Int_t nStart, Int_t nStop, Bool_t copyCache) :
00888 RooAbsData(name,title,varSubset),
00889 _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10)
00890 {
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901 _dstore = new RooTreeDataStore(name,title,*h->_dstore,_vars,cutVar,cutRange,nStart,nStop,copyCache) ;
00902
00903 initialize(0,kFALSE) ;
00904
00905 ((RooTreeDataStore*)_dstore)->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
00906
00907
00908 Int_t i ;
00909 for (i=0 ; i<_arrSize ; i++) {
00910 _wgt[i] = h->_wgt[i] ;
00911 _errLo[i] = h->_errLo[i] ;
00912 _errHi[i] = h->_errHi[i] ;
00913 _sumw2[i] = h->_sumw2[i] ;
00914 _binv[i] = h->_binv[i] ;
00915 }
00916
00917
00918 appendToDir(this,kTRUE) ;
00919 }
00920
00921
00922
00923 RooAbsData* RooDataHist::cacheClone(const RooAbsArg* newCacheOwner, const RooArgSet* newCacheVars, const char* newName)
00924 {
00925
00926 checkInit() ;
00927
00928 RooDataHist* dhist = new RooDataHist(newName?newName:GetName(),GetTitle(),this,*get(),0,0,0,2000000000,kTRUE) ;
00929
00930 RooArgSet* selCacheVars = (RooArgSet*) newCacheVars->selectCommon(dhist->_cachedVars) ;
00931 dhist->attachCache(newCacheOwner, *selCacheVars) ;
00932 delete selCacheVars ;
00933
00934 return dhist ;
00935 }
00936
00937
00938
00939
00940 RooAbsData* RooDataHist::reduceEng(const RooArgSet& varSubset, const RooFormulaVar* cutVar, const char* cutRange,
00941 Int_t nStart, Int_t nStop, Bool_t )
00942 {
00943
00944
00945 checkInit() ;
00946 RooArgSet* myVarSubset = (RooArgSet*) _vars.selectCommon(varSubset) ;
00947 RooDataHist *rdh = new RooDataHist(GetName(), GetTitle(), *myVarSubset) ;
00948 delete myVarSubset ;
00949
00950 RooFormulaVar* cloneVar = 0;
00951 RooArgSet* tmp(0) ;
00952 if (cutVar) {
00953
00954 tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
00955 if (!tmp) {
00956 coutE(DataHandling) << "RooDataHist::reduceEng(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
00957 return 0 ;
00958 }
00959 cloneVar = (RooFormulaVar*) tmp->find(cutVar->GetName()) ;
00960 cloneVar->attachDataSet(*this) ;
00961 }
00962
00963 Int_t i ;
00964 Double_t lo,hi ;
00965 Int_t nevt = nStop < numEntries() ? nStop : numEntries() ;
00966 TIterator* vIter = get()->createIterator() ;
00967 for (i=nStart ; i<nevt ; i++) {
00968 const RooArgSet* row = get(i) ;
00969
00970 Bool_t doSelect(kTRUE) ;
00971 if (cutRange) {
00972 RooAbsArg* arg ;
00973 vIter->Reset() ;
00974 while((arg=(RooAbsArg*)vIter->Next())) {
00975 if (!arg->inRange(cutRange)) {
00976 doSelect = kFALSE ;
00977 break ;
00978 }
00979 }
00980 }
00981 if (!doSelect) continue ;
00982
00983 if (!cloneVar || cloneVar->getVal()) {
00984 weightError(lo,hi,SumW2) ;
00985 rdh->add(*row,weight(),lo*lo) ;
00986 }
00987 }
00988 delete vIter ;
00989
00990 if (cloneVar) {
00991 delete tmp ;
00992 }
00993
00994 return rdh ;
00995 }
00996
00997
00998
00999
01000 RooDataHist::~RooDataHist()
01001 {
01002
01003
01004 if (_wgt) delete[] _wgt ;
01005 if (_errLo) delete[] _errLo ;
01006 if (_errHi) delete[] _errHi ;
01007 if (_sumw2) delete[] _sumw2 ;
01008 if (_binv) delete[] _binv ;
01009 if (_realIter) delete _realIter ;
01010 if (_binValid) delete[] _binValid ;
01011 list<const RooAbsBinning*>::iterator iter = _lvbins.begin() ;
01012 while(iter!=_lvbins.end()) {
01013 delete *iter ;
01014 iter++ ;
01015 }
01016
01017 removeFromDir(this) ;
01018 }
01019
01020
01021
01022
01023
01024 Int_t RooDataHist::getIndex(const RooArgSet& coord)
01025 {
01026 checkInit() ;
01027 _vars = coord ;
01028 return calcTreeIndex() ;
01029 }
01030
01031
01032
01033
01034
01035 Int_t RooDataHist::calcTreeIndex() const
01036 {
01037
01038
01039
01040 Int_t masterIdx(0), i(0) ;
01041 list<RooAbsLValue*>::const_iterator iter = _lvvars.begin() ;
01042 list<const RooAbsBinning*>::const_iterator biter = _lvbins.begin() ;
01043 for (;iter!=_lvvars.end() ; ++iter) {
01044 const RooAbsBinning* binning = (*biter) ;
01045 masterIdx += _idxMult[i++]*(*iter)->getBin(binning) ;
01046 biter++ ;
01047 }
01048 return masterIdx ;
01049 }
01050
01051
01052
01053
01054 void RooDataHist::dump2()
01055 {
01056
01057 cout << "_arrSize = " << _arrSize << endl ;
01058 for (Int_t i=0 ; i<_arrSize ; i++) {
01059 cout << "wgt[" << i << "] = " << _wgt[i] << "sumw2[" << i << "] = " << _sumw2[i] << " vol[" << i << "] = " << _binv[i] << endl ;
01060 }
01061 }
01062
01063
01064
01065
01066 RooPlot *RooDataHist::plotOn(RooPlot *frame, PlotOpt o) const
01067 {
01068
01069
01070
01071
01072 checkInit() ;
01073 if (o.bins) return RooAbsData::plotOn(frame,o) ;
01074
01075 if(0 == frame) {
01076 coutE(InputArguments) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
01077 return 0;
01078 }
01079 RooAbsRealLValue *var= (RooAbsRealLValue*) frame->getPlotVar();
01080 if(0 == var) {
01081 coutE(InputArguments) << ClassName() << "::" << GetName()
01082 << ":plotOn: frame does not specify a plot variable" << endl;
01083 return 0;
01084 }
01085
01086 RooRealVar* dataVar = (RooRealVar*) _vars.find(var->GetName()) ;
01087 if (!dataVar) {
01088 coutE(InputArguments) << ClassName() << "::" << GetName()
01089 << ":plotOn: dataset doesn't contain plot frame variable" << endl;
01090 return 0;
01091 }
01092
01093 o.bins = &dataVar->getBinning() ;
01094 o.correctForBinWidth = kFALSE ;
01095 return RooAbsData::plotOn(frame,o) ;
01096 }
01097
01098
01099
01100
01101
01102 Double_t RooDataHist::weight(const RooArgSet& bin, Int_t intOrder, Bool_t correctForBinSize, Bool_t cdfBoundaries)
01103 {
01104
01105
01106
01107
01108
01109
01110
01111 checkInit() ;
01112
01113
01114 if (intOrder<0) {
01115 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
01116 return 0 ;
01117 }
01118
01119
01120 if (intOrder==0) {
01121 _vars.assignValueOnly(bin) ;
01122 Int_t idx = calcTreeIndex() ;
01123
01124 if (correctForBinSize) {
01125
01126 return _wgt[idx] / _binv[idx] ;
01127 } else {
01128 return _wgt[idx] ;
01129 }
01130 }
01131
01132
01133 _vars.assignValueOnly(bin) ;
01134
01135 Double_t wInt(0) ;
01136 if (_realVars.getSize()==1) {
01137
01138
01139 _realIter->Reset() ;
01140 RooRealVar* real=(RooRealVar*)_realIter->Next() ;
01141 const RooAbsBinning* binning = real->getBinningPtr(0) ;
01142 wInt = interpolateDim(*real,binning,((RooAbsReal*)bin.find(real->GetName()))->getVal(), intOrder, correctForBinSize, cdfBoundaries) ;
01143
01144 } else if (_realVars.getSize()==2) {
01145
01146
01147 _realIter->Reset() ;
01148 RooRealVar* realX=(RooRealVar*)_realIter->Next() ;
01149 RooRealVar* realY=(RooRealVar*)_realIter->Next() ;
01150 Double_t xval = ((RooAbsReal*)bin.find(realX->GetName()))->getVal() ;
01151 Double_t yval = ((RooAbsReal*)bin.find(realY->GetName()))->getVal() ;
01152
01153 Int_t ybinC = realY->getBin() ;
01154 Int_t ybinLo = ybinC-intOrder/2 - ((yval<realY->getBinning().binCenter(ybinC))?1:0) ;
01155 Int_t ybinM = realY->numBins() ;
01156
01157 Int_t i ;
01158 Double_t yarr[10] ;
01159 Double_t xarr[10] ;
01160 const RooAbsBinning* binning = realX->getBinningPtr(0) ;
01161 for (i=ybinLo ; i<=intOrder+ybinLo ; i++) {
01162 Int_t ibin ;
01163 if (i>=0 && i<ybinM) {
01164
01165 ibin = i ;
01166 realY->setBin(ibin) ;
01167 xarr[i-ybinLo] = realY->getVal() ;
01168 } else if (i>=ybinM) {
01169
01170 ibin = 2*ybinM-i-1 ;
01171 realY->setBin(ibin) ;
01172 xarr[i-ybinLo] = 2*realY->getMax()-realY->getVal() ;
01173 } else {
01174
01175 ibin = -i ;
01176 realY->setBin(ibin) ;
01177 xarr[i-ybinLo] = 2*realY->getMin()-realY->getVal() ;
01178 }
01179 yarr[i-ybinLo] = interpolateDim(*realX,binning,xval,intOrder,correctForBinSize,kFALSE) ;
01180 }
01181
01182 if (gDebug>7) {
01183 cout << "RooDataHist interpolating data is" << endl ;
01184 cout << "xarr = " ;
01185 for (int q=0; q<=intOrder ; q++) cout << xarr[q] << " " ;
01186 cout << " yarr = " ;
01187 for (int q=0; q<=intOrder ; q++) cout << yarr[q] << " " ;
01188 cout << endl ;
01189 }
01190 wInt = RooMath::interpolate(xarr,yarr,intOrder+1,yval) ;
01191
01192 } else {
01193
01194
01195 coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") interpolation in "
01196 << _realVars.getSize() << " dimensions not yet implemented" << endl ;
01197 return weight(bin,0) ;
01198
01199 }
01200
01201
01202
01203
01204
01205
01206
01207 return wInt ;
01208 }
01209
01210
01211
01212
01213
01214 void RooDataHist::weightError(Double_t& lo, Double_t& hi, ErrorType etype) const
01215 {
01216
01217 checkInit() ;
01218
01219 switch (etype) {
01220
01221 case Auto:
01222 throw string(Form("RooDataHist::weightError(%s) weight type Auto not allowed here",GetName())) ;
01223 break ;
01224
01225 case Poisson:
01226 if (_curWgtErrLo>=0) {
01227
01228 lo = _curWgtErrLo ;
01229 hi = _curWgtErrHi ;
01230 return ;
01231 }
01232
01233
01234 Double_t ym,yp ;
01235 RooHistError::instance().getPoissonInterval(Int_t(weight()+0.5),ym,yp,1) ;
01236 _curWgtErrLo = weight()-ym ;
01237 _curWgtErrHi = yp-weight() ;
01238 _errLo[_curIndex] = _curWgtErrLo ;
01239 _errHi[_curIndex] = _curWgtErrHi ;
01240 lo = _curWgtErrLo ;
01241 hi = _curWgtErrHi ;
01242 return ;
01243
01244 case SumW2:
01245 lo = sqrt(_curSumW2) ;
01246 hi = sqrt(_curSumW2) ;
01247 return ;
01248
01249 case None:
01250 lo = 0 ;
01251 hi = 0 ;
01252 return ;
01253 }
01254 }
01255
01256
01257
01258
01259
01260 Double_t RooDataHist::interpolateDim(RooRealVar& dim, const RooAbsBinning* binning, Double_t xval, Int_t intOrder, Bool_t correctForBinSize, Bool_t cdfBoundaries)
01261 {
01262
01263
01264
01265
01266 Int_t fbinC = dim.getBin(*binning) ;
01267 Int_t fbinLo = fbinC-intOrder/2 - ((xval<binning->binCenter(fbinC))?1:0) ;
01268 Int_t fbinM = dim.numBins(*binning) ;
01269
01270
01271 Int_t i ;
01272 Double_t yarr[10] ;
01273 Double_t xarr[10] ;
01274 for (i=fbinLo ; i<=intOrder+fbinLo ; i++) {
01275 Int_t ibin ;
01276 if (i>=0 && i<fbinM) {
01277
01278 ibin = i ;
01279 dim.setBinFast(ibin,*binning) ;
01280
01281 xarr[i-fbinLo] = dim.getVal() ;
01282 Int_t idx = calcTreeIndex() ;
01283 yarr[i-fbinLo] = _wgt[idx] ;
01284 if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
01285 } else if (i>=fbinM) {
01286
01287 ibin = 2*fbinM-i-1 ;
01288 dim.setBinFast(ibin,*binning) ;
01289
01290 if (cdfBoundaries) {
01291 xarr[i-fbinLo] = dim.getMax()+1e-10*(i-fbinM+1) ;
01292 yarr[i-fbinLo] = 1.0 ;
01293 } else {
01294 Int_t idx = calcTreeIndex() ;
01295 xarr[i-fbinLo] = 2*dim.getMax()-dim.getVal() ;
01296 yarr[i-fbinLo] = _wgt[idx] ;
01297 if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
01298 }
01299 } else {
01300
01301 ibin = -i - 1 ;
01302 dim.setBinFast(ibin,*binning) ;
01303
01304 if (cdfBoundaries) {
01305 xarr[i-fbinLo] = dim.getMin()-ibin*(1e-10) ; ;
01306 yarr[i-fbinLo] = 0.0 ;
01307 } else {
01308 Int_t idx = calcTreeIndex() ;
01309 xarr[i-fbinLo] = 2*dim.getMin()-dim.getVal() ;
01310 yarr[i-fbinLo] = _wgt[idx] ;
01311 if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
01312 }
01313 }
01314
01315 }
01316
01317
01318
01319 dim.setBinFast(fbinC,*binning) ;
01320 Double_t ret = RooMath::interpolate(xarr,yarr,intOrder+1,xval) ;
01321 return ret ;
01322 }
01323
01324
01325
01326
01327
01328 void RooDataHist::add(const RooArgSet& row, Double_t wgt, Double_t sumw2)
01329 {
01330
01331
01332
01333
01334 checkInit() ;
01335
01336
01337
01338
01339 _vars = row ;
01340 Int_t idx = calcTreeIndex() ;
01341 _wgt[idx] += wgt ;
01342 _sumw2[idx] += (sumw2>0?sumw2:wgt*wgt) ;
01343 _errLo[idx] = -1 ;
01344 _errHi[idx] = -1 ;
01345 }
01346
01347
01348
01349
01350 void RooDataHist::set(const RooArgSet& row, Double_t wgt, Double_t wgtErrLo, Double_t wgtErrHi)
01351 {
01352
01353
01354
01355
01356 checkInit() ;
01357
01358 _vars = row ;
01359 Int_t idx = calcTreeIndex() ;
01360 _wgt[idx] = wgt ;
01361 _errLo[idx] = wgtErrLo ;
01362 _errHi[idx] = wgtErrHi ;
01363 }
01364
01365
01366
01367
01368 void RooDataHist::set(Double_t wgt, Double_t wgtErr)
01369 {
01370
01371
01372
01373
01374 checkInit() ;
01375
01376 if (_curIndex<0) {
01377 _curIndex = calcTreeIndex() ;
01378 }
01379
01380 _wgt[_curIndex] = wgt ;
01381 _errLo[_curIndex] = wgtErr ;
01382 _errHi[_curIndex] = wgtErr ;
01383 _sumw2[_curIndex] = wgtErr*wgtErr ;
01384 }
01385
01386
01387
01388
01389 void RooDataHist::set(const RooArgSet& row, Double_t wgt, Double_t wgtErr)
01390 {
01391
01392
01393
01394
01395 checkInit() ;
01396
01397 _vars = row ;
01398 Int_t idx = calcTreeIndex() ;
01399 _wgt[idx] = wgt ;
01400 _errLo[idx] = wgtErr ;
01401 _errHi[idx] = wgtErr ;
01402 _sumw2[idx] = wgtErr*wgtErr ;
01403 }
01404
01405
01406
01407
01408 void RooDataHist::add(const RooAbsData& dset, const char* cut, Double_t wgt)
01409 {
01410
01411
01412
01413
01414 RooFormulaVar cutVar("select",cut,*dset.get()) ;
01415 add(dset,&cutVar,wgt) ;
01416 }
01417
01418
01419
01420
01421 void RooDataHist::add(const RooAbsData& dset, const RooFormulaVar* cutVar, Double_t wgt)
01422 {
01423
01424
01425
01426 checkInit() ;
01427
01428 RooFormulaVar* cloneVar = 0;
01429 RooArgSet* tmp(0) ;
01430 if (cutVar) {
01431
01432 tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
01433 if (!tmp) {
01434 coutE(DataHandling) << "RooDataHist::add(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
01435 return ;
01436 }
01437
01438 cloneVar = (RooFormulaVar*) tmp->find(cutVar->GetName()) ;
01439 cloneVar->attachDataSet(dset) ;
01440 }
01441
01442
01443 Int_t i ;
01444 for (i=0 ; i<dset.numEntries() ; i++) {
01445 const RooArgSet* row = dset.get(i) ;
01446 if (!cloneVar || cloneVar->getVal()) {
01447 add(*row,wgt*dset.weight()) ;
01448 }
01449 }
01450
01451 if (cloneVar) {
01452 delete tmp ;
01453 }
01454 }
01455
01456
01457
01458
01459 Double_t RooDataHist::sum(Bool_t correctForBinSize) const
01460 {
01461
01462
01463
01464
01465
01466
01467
01468 checkInit() ;
01469
01470 Int_t i ;
01471 Double_t total(0) ;
01472 for (i=0 ; i<_arrSize ; i++) {
01473
01474 Double_t theBinVolume = correctForBinSize ? _binv[i] : 1.0 ;
01475
01476 total += _wgt[i]*theBinVolume ;
01477 }
01478
01479 return total ;
01480 }
01481
01482
01483
01484
01485 Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bool_t correctForBinSize)
01486 {
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497 checkInit() ;
01498
01499 RooArgSet varSave ;
01500 varSave.addClone(_vars) ;
01501
01502 RooArgSet* sliceOnlySet = new RooArgSet(sliceSet) ;
01503 sliceOnlySet->remove(sumSet,kTRUE,kTRUE) ;
01504
01505 _vars = *sliceOnlySet ;
01506 calculatePartialBinVolume(*sliceOnlySet) ;
01507 delete sliceOnlySet ;
01508
01509 TIterator* ssIter = sumSet.createIterator() ;
01510
01511
01512 RooAbsArg* arg ;
01513 Bool_t* mask = new Bool_t[_vars.getSize()] ;
01514 Int_t* refBin = new Int_t[_vars.getSize()] ;
01515
01516 Int_t i(0) ;
01517 _iterator->Reset() ;
01518 while((arg=(RooAbsArg*)_iterator->Next())) {
01519 if (sumSet.find(arg->GetName())) {
01520 mask[i] = kFALSE ;
01521 } else {
01522 mask[i] = kTRUE ;
01523
01524 refBin[i] = (dynamic_cast<RooAbsLValue*>(arg))->getBin() ;
01525 }
01526 i++ ;
01527 }
01528
01529
01530 Double_t total(0) ;
01531 Int_t ibin ;
01532 for (ibin=0 ; ibin<_arrSize ; ibin++) {
01533
01534 Int_t idx(0), tmp(ibin), ivar(0) ;
01535 Bool_t skip(kFALSE) ;
01536
01537
01538 _iterator->Reset() ;
01539
01540 while((!skip && (arg=(RooAbsArg*)_iterator->Next()))) {
01541 idx = tmp / _idxMult[ivar] ;
01542 tmp -= idx*_idxMult[ivar] ;
01543 if (mask[ivar] && idx!=refBin[ivar]) skip=kTRUE ;
01544 ivar++ ;
01545 }
01546
01547 if (!skip) {
01548 Double_t theBinVolume = correctForBinSize ? (*_pbinv)[ibin] : 1.0 ;
01549
01550 total += _wgt[ibin]/theBinVolume ;
01551 }
01552 }
01553 delete ssIter ;
01554
01555 delete[] mask ;
01556 delete[] refBin ;
01557
01558 _vars = varSave ;
01559
01560 return total ;
01561 }
01562
01563
01564
01565
01566 void RooDataHist::calculatePartialBinVolume(const RooArgSet& dimSet) const
01567 {
01568
01569
01570
01571
01572 vector<Double_t> *pbinv = _pbinvCacheMgr.getObj(&dimSet) ;
01573 if (pbinv) {
01574 _pbinv = pbinv ;
01575 return ;
01576 }
01577
01578 pbinv = new vector<Double_t>(_arrSize) ;
01579
01580
01581 Bool_t* selDim = new Bool_t[_vars.getSize()] ;
01582 _iterator->Reset() ;
01583 RooAbsArg* v ;
01584 Int_t i(0) ;
01585 while((v=(RooAbsArg*)_iterator->Next())) {
01586 selDim[i++] = dimSet.find(v->GetName()) ? kTRUE : kFALSE ;
01587 }
01588
01589
01590 Int_t ibin ;
01591 for (ibin=0 ; ibin<_arrSize ; ibin++) {
01592 _iterator->Reset() ;
01593 RooAbsLValue* arg ;
01594 Int_t j(0), idx(0), tmp(ibin) ;
01595 Double_t theBinVolume(1) ;
01596 while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
01597 idx = tmp / _idxMult[j] ;
01598 tmp -= idx*_idxMult[j++] ;
01599 if (selDim[j-1]) {
01600 RooAbsLValue* arglv = dynamic_cast<RooAbsLValue*>(arg) ;
01601 theBinVolume *= arglv->getBinWidth(idx) ;
01602 }
01603 }
01604 (*pbinv)[ibin] = theBinVolume ;
01605 }
01606
01607 delete[] selDim ;
01608
01609
01610 _pbinvCacheMgr.setObj(&dimSet,pbinv) ;
01611
01612
01613 _pbinv = pbinv ;
01614 }
01615
01616
01617
01618
01619 Int_t RooDataHist::numEntries() const
01620 {
01621
01622
01623 return RooAbsData::numEntries() ;
01624 }
01625
01626
01627
01628
01629 Double_t RooDataHist::sumEntries(const char* cutSpec, const char* cutRange) const
01630 {
01631
01632
01633
01634 checkInit() ;
01635 if (cutSpec==0 && cutRange==0) {
01636 Int_t i ;
01637 Double_t n(0) ;
01638 for (i=0 ; i<_arrSize ; i++) {
01639 if (!_binValid || _binValid[i]) {
01640 n+= _wgt[i] ;
01641 }
01642 }
01643 return n ;
01644 } else {
01645
01646
01647 RooFormula* select = 0 ;
01648 if (cutSpec) {
01649 select = new RooFormula("select",cutSpec,*get()) ;
01650 }
01651
01652
01653 Double_t sumw(0) ;
01654 Int_t i ;
01655 for (i=0 ; i<numEntries() ; i++) {
01656 get(i) ;
01657 if (select && select->eval()==0.) continue ;
01658 if (cutRange && !_vars.allInRange(cutRange)) continue ;
01659
01660 if (!_binValid || _binValid[i]) {
01661 sumw += weight() ;
01662 }
01663 }
01664
01665 if (select) delete select ;
01666
01667 return sumw ;
01668 }
01669 }
01670
01671
01672
01673
01674 void RooDataHist::reset()
01675 {
01676
01677
01678
01679
01680
01681 Int_t i ;
01682 for (i=0 ; i<_arrSize ; i++) {
01683 _wgt[i] = 0. ;
01684 _errLo[i] = -1 ;
01685 _errHi[i] = -1 ;
01686 }
01687 _curWeight = 0 ;
01688 _curWgtErrLo = -1 ;
01689 _curWgtErrHi = -1 ;
01690 _curVolume = 1 ;
01691
01692 }
01693
01694
01695
01696
01697 const RooArgSet* RooDataHist::get(Int_t masterIdx) const
01698 {
01699
01700
01701 checkInit() ;
01702 _curWeight = _wgt[masterIdx] ;
01703 _curWgtErrLo = _errLo[masterIdx] ;
01704 _curWgtErrHi = _errHi[masterIdx] ;
01705 _curSumW2 = _sumw2[masterIdx] ;
01706 _curVolume = _binv[masterIdx] ;
01707 _curIndex = masterIdx ;
01708 return RooAbsData::get(masterIdx) ;
01709 }
01710
01711
01712
01713
01714 const RooArgSet* RooDataHist::get(const RooArgSet& coord) const
01715 {
01716
01717
01718
01719 ((RooDataHist*)this)->_vars = coord ;
01720 return get(calcTreeIndex()) ;
01721 }
01722
01723
01724
01725
01726 Double_t RooDataHist::binVolume(const RooArgSet& coord)
01727 {
01728
01729 checkInit() ;
01730 ((RooDataHist*)this)->_vars = coord ;
01731 return _binv[calcTreeIndex()] ;
01732 }
01733
01734
01735
01736 void RooDataHist::setAllWeights(Double_t value)
01737 {
01738
01739 for (Int_t i=0 ; i<_arrSize ; i++) {
01740 _wgt[i] = value ;
01741 }
01742 }
01743
01744
01745
01746
01747 TIterator* RooDataHist::sliceIterator(RooAbsArg& sliceArg, const RooArgSet& otherArgs)
01748 {
01749
01750
01751
01752
01753 _vars = otherArgs ;
01754 _curIndex = calcTreeIndex() ;
01755
01756 RooAbsArg* intArg = _vars.find(sliceArg.GetName()) ;
01757 if (!intArg) {
01758 coutE(InputArguments) << "RooDataHist::sliceIterator() variable " << sliceArg.GetName() << " is not part of this RooDataHist" << endl ;
01759 return 0 ;
01760 }
01761 return new RooDataHistSliceIter(*this,*intArg) ;
01762 }
01763
01764
01765
01766 void RooDataHist::SetName(const char *name)
01767 {
01768
01769
01770 if (_dir) _dir->GetList()->Remove(this);
01771 TNamed::SetName(name) ;
01772 if (_dir) _dir->GetList()->Add(this);
01773 }
01774
01775
01776
01777 void RooDataHist::SetNameTitle(const char *name, const char* title)
01778 {
01779
01780
01781 if (_dir) _dir->GetList()->Remove(this);
01782 TNamed::SetNameTitle(name,title) ;
01783 if (_dir) _dir->GetList()->Add(this);
01784 }
01785
01786
01787
01788 void RooDataHist::printValue(ostream& os) const
01789 {
01790
01791 os << numEntries() << " bins (" << sumEntries() << " weights)" ;
01792 }
01793
01794
01795
01796
01797
01798 void RooDataHist::printArgs(ostream& os) const
01799 {
01800
01801
01802 os << "[" ;
01803 _iterator->Reset() ;
01804 RooAbsArg* arg ;
01805 Bool_t first(kTRUE) ;
01806 while((arg=(RooAbsArg*)_iterator->Next())) {
01807 if (first) {
01808 first=kFALSE ;
01809 } else {
01810 os << "," ;
01811 }
01812 os << arg->GetName() ;
01813 }
01814 os << "]" ;
01815 }
01816
01817
01818
01819
01820 void RooDataHist::cacheValidEntries()
01821 {
01822
01823
01824 checkInit() ;
01825
01826 if (!_binValid) {
01827 _binValid = new Bool_t[_arrSize] ;
01828 }
01829 TIterator* iter = _vars.createIterator() ;
01830 RooAbsArg* arg ;
01831 for (Int_t i=0 ; i<_arrSize ; i++) {
01832 get(i) ;
01833 _binValid[i] = kTRUE ;
01834 iter->Reset() ;
01835 while((arg=(RooAbsArg*)iter->Next())) {
01836
01837 _binValid[i] &= arg->inRange(0) ;
01838 }
01839 }
01840 delete iter ;
01841
01842 }
01843
01844
01845
01846 Bool_t RooDataHist::valid() const
01847 {
01848
01849
01850
01851
01852 if (_binValid) {
01853 return _binValid[_curIndex] ;
01854 }
01855
01856 return kTRUE ;
01857 }
01858
01859
01860
01861
01862 Bool_t RooDataHist::isNonPoissonWeighted() const
01863 {
01864
01865
01866 for (int i=0 ; i<numEntries() ; i++) {
01867 if (fabs(_wgt[i]-Int_t(_wgt[i]))>1e-10) return kTRUE ;
01868 }
01869 return kFALSE ;
01870 }
01871
01872
01873
01874
01875
01876 void RooDataHist::printMultiline(ostream& os, Int_t content, Bool_t verbose, TString indent) const
01877 {
01878
01879
01880 RooAbsData::printMultiline(os,content,verbose,indent) ;
01881
01882 os << indent << "Binned Dataset " << GetName() << " (" << GetTitle() << ")" << endl ;
01883 os << indent << " Contains " << numEntries() << " bins with a total weight of " << sumEntries() << endl;
01884
01885 if (!verbose) {
01886 os << indent << " Observables " << _vars << endl ;
01887 } else {
01888 os << indent << " Observables: " ;
01889 _vars.printStream(os,kName|kValue|kExtras|kTitle,kVerbose,indent+" ") ;
01890 }
01891
01892 if(verbose) {
01893 if (_cachedVars.getSize()>0) {
01894 os << indent << " Caches " << _cachedVars << endl ;
01895 }
01896 }
01897 }
01898
01899
01900
01901
01902 void RooDataHist::Streamer(TBuffer &R__b)
01903 {
01904
01905
01906 if (R__b.IsReading()) {
01907
01908 UInt_t R__s, R__c;
01909 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
01910
01911 if (R__v>2) {
01912
01913 R__b.ReadClassBuffer(RooDataHist::Class(),this,R__v,R__s,R__c);
01914 initialize(0,kFALSE) ;
01915
01916 } else {
01917
01918
01919
01920
01921
01922
01923
01924 UInt_t R__s1, R__c1;
01925 Version_t R__v1 = R__b.ReadVersion(&R__s1, &R__c1); if (R__v1) { }
01926
01927 RooAbsData::Streamer(R__b);
01928 TTree* X_tree(0) ; R__b >> X_tree;
01929 RooArgSet X_truth ; X_truth.Streamer(R__b);
01930 TString X_blindString ; X_blindString.Streamer(R__b);
01931 R__b.CheckByteCount(R__s1, R__c1, RooTreeData::Class());
01932
01933
01934
01935 _dstore = new RooTreeDataStore(X_tree,_vars) ;
01936 _dstore->SetName(GetName()) ;
01937 _dstore->SetTitle(GetTitle()) ;
01938 _dstore->checkInit() ;
01939
01940 RooDirItem::Streamer(R__b);
01941 R__b >> _arrSize;
01942 delete [] _wgt;
01943 _wgt = new Double_t[_arrSize];
01944 R__b.ReadFastArray(_wgt,_arrSize);
01945 delete [] _errLo;
01946 _errLo = new Double_t[_arrSize];
01947 R__b.ReadFastArray(_errLo,_arrSize);
01948 delete [] _errHi;
01949 _errHi = new Double_t[_arrSize];
01950 R__b.ReadFastArray(_errHi,_arrSize);
01951 delete [] _sumw2;
01952 _sumw2 = new Double_t[_arrSize];
01953 R__b.ReadFastArray(_sumw2,_arrSize);
01954 delete [] _binv;
01955 _binv = new Double_t[_arrSize];
01956 R__b.ReadFastArray(_binv,_arrSize);
01957 _realVars.Streamer(R__b);
01958 R__b >> _curWeight;
01959 R__b >> _curWgtErrLo;
01960 R__b >> _curWgtErrHi;
01961 R__b >> _curSumW2;
01962 R__b >> _curVolume;
01963 R__b >> _curIndex;
01964 R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
01965
01966 }
01967
01968 } else {
01969
01970 R__b.WriteClassBuffer(RooDataHist::Class(),this);
01971 }
01972 }
01973
01974