119 fFoamType(kSeparate),
126 fFillFoamWithOrigWeights(
kFALSE),
127 fDTSeparation(kFoam),
158 fFoamType(kSeparate),
165 fFillFoamWithOrigWeights(
kFALSE),
166 fDTSeparation(kFoam),
173 if(strlen(
name) > 128)
174 Log() << kFATAL <<
"Name too long " <<
name.Data() <<
Endl;
186 delete fVariableNames;
188 if (fDistr)
delete fDistr;
189 if (fPseRan)
delete fPseRan;
190 if (fXmin) {
delete [] fXmin; fXmin=0; }
191 if (fXmax) {
delete [] fXmax; fXmax=0; }
195 for(
Int_t i=0; i<fNCells; i++)
delete fCells[i];
225 , fFoamType(kSeparate)
232 , fFillFoamWithOrigWeights(
kFALSE)
233 , fDTSeparation(kFoam)
240 Log() << kFATAL <<
"COPY CONSTRUCTOR NOT IMPLEMENTED" <<
Endl;
253 Log() << kFATAL <<
"<SetDim>: Dimension is zero or negative!" <<
Endl;
256 if (fXmin)
delete [] fXmin;
257 if (fXmax)
delete [] fXmax;
267 if (idim<0 || idim>=GetTotDim())
268 Log() << kFATAL <<
"<SetXmin>: Dimension out of bounds!" <<
Endl;
278 if (idim<0 || idim>=GetTotDim())
279 Log() << kFATAL <<
"<SetXmax>: Dimension out of bounds!" <<
Endl;
297 if(fPseRan==0)
Log() << kFATAL <<
"Random number generator not set" <<
Endl;
298 if(fDistr==0)
Log() << kFATAL <<
"Distribution function not set" <<
Endl;
299 if(fDim==0)
Log() << kFATAL <<
"Zero dimension not allowed" <<
Endl;
306 if(fRvec==0)
Log() << kFATAL <<
"Cannot initialize buffer fRvec" <<
Endl;
310 if(fAlpha==0)
Log() << kFATAL <<
"Cannot initialize buffer fAlpha" <<
Endl;
315 fInhiDiv =
new Int_t[fDim];
316 for(
Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
320 fMaskDiv =
new Int_t[fDim];
321 for(
Int_t i=0; i<fDim; i++) fMaskDiv[i]=1;
325 for(
Int_t i=0; i<fDim; i++){
327 hname = fName+
TString(
"_HistEdge_");
329 htitle =
TString(
"Edge Histogram No. ");
331 (*fHistEdg)[i] =
new TH1D(hname.
Data(),htitle.
Data(),fNBin,0.0, 1.0);
332 ((
TH1D*)(*fHistEdg)[i])->Sumw2();
360 for(
Int_t i=0; i<fNCells; i++)
delete fCells[i];
366 Log() << kFATAL <<
"not enough memory to create " << fNCells
369 for(
Int_t i=0; i<fNCells; i++){
371 fCells[i]->SetSerial(i);
380 for(
Long_t iCell=0; iCell<=fLastCe; iCell++){
381 Explore( fCells[iCell] );
392 if (fLastCe==fNCells){
393 Log() << kFATAL <<
"Too many cells" <<
Endl;
397 cell = fCells[fLastCe];
399 cell->
Fill(status, parent, 0, 0);
452 cell->
GetHcub(cellPosi,cellSize);
462 for (
Int_t idim = 0; idim < fDim; ++idim)
463 vol_scale *= fXmax[idim] - fXmin[idim];
469 toteventsOld = GetCellElement(cell, 0);
480 for (i=0;i<fDim;i++) ((
TH1D *)(*fHistEdg)[i])->Reset();
484 for (iev=0;iev<fNSampl;iev++){
487 if (fDim>0)
for (j=0; j<fDim; j++) xRand[j]= cellPosi[j] +fAlpha[j]*(cellSize[j]);
489 wt = dx*
Eval(xRand, event_density);
490 totevents += event_density;
494 for (k=0; k<fDim; k++) {
496 ((
TH1D *)(*fHistEdg)[nProj])->Fill(xproj,wt);
504 if (ceSum[3]>wt) ceSum[3]=wt;
505 if (ceSum[4]<wt) ceSum[4]=wt;
507 if (ceSum[1]>0) nevEff = ceSum[0]*ceSum[0]/ceSum[1];
509 if ( nevEff >= fNBin*fEvPerBin)
break;
513 if (fNSampl>0) totevents /= fNSampl;
517 if (cell==fCells[0] && ceSum[0]<=0.0){
519 Log() << kFATAL <<
"No events were found during exploration of "
520 <<
"root cell. Please check PDEFoam parameters nSampl "
521 <<
"and VolFrac." <<
Endl;
523 Log() << kWARNING <<
"Negative number of events found during "
524 <<
"exploration of root cell" <<
Endl;
529 for (k=0; k<fDim;k++){
531 if ( fInhiDiv[k]==1) fMaskDiv[k] =0;
538 Double_t intTrue = ceSum[0]/(nevMC+0.000001);
541 if (kBest == -1) Varedu(ceSum,kBest,xBest,yBest);
542 intDriv =
sqrt(ceSum[1]/nevMC) -intTrue;
549 SetCellElement(cell, 0, totevents);
553 for (parent = cell->
GetPare(); parent!=0; parent = parent->
GetPare()){
556 parent->
SetIntg( parIntg +intTrue -intOld );
557 parent->
SetDriv( parDriv +intDriv -driOld );
558 SetCellElement( parent, 0, GetCellElement(parent, 0) + totevents - toteventsOld);
583 for(
Int_t kProj=0; kProj<fDim; kProj++) {
584 if( fMaskDiv[kProj]) {
591 for(
Int_t jLo=1; jLo<=fNBin; jLo++) {
593 for(
Int_t jUp=jLo; jUp<=fNBin;jUp++) {
594 aswIn += ((
TH1D *)(*fHistEdg)[kProj])->GetBinContent(jUp);
595 asswIn += Sqr(((
TH1D *)(*fHistEdg)[kProj])->GetBinError( jUp));
601 else sswIn =
sqrt(asswIn) /
sqrt(nent*(xUp-xLo)) *(xUp-xLo);
604 else sswOut=
sqrt(sswAll-asswIn)/
sqrt(nent*(1.0-xUp+xLo)) *(1.0-xUp+xLo);
605 if( (sswIn+sswOut) < sswtBest) {
606 sswtBest = sswIn+sswOut;
623 if(iLo == 0) xBest=yBest;
624 if(iUp == fNBin) yBest=xBest;
629 if( (kBest >= fDim) || (kBest<0) )
630 Log() << kFATAL <<
"Something wrong with kBest" <<
Endl;
640 fPseRan->RndmArray(fDim,fRvec);
641 for(
Int_t k=0; k<fDim; k++) fAlpha[k] = fRvec[k];
661 for(i=0; i<=fLastCe; i++) {
662 if( fCells[i]->GetStat() == 1 ) {
664 driv = fCells[i]->GetDriv();
675 if (GetMaxDepth() > 0)
676 bCutMaxDepth = fCells[i]->GetDepth() < GetMaxDepth();
680 bCutNmin = GetCellElement(fCells[i], 0) > GetNmin();
683 if(driv > drivMax && bCutNmin && bCutMaxDepth) {
692 Log() << kVERBOSE <<
"Warning: No cell with more than "
693 << GetNmin() <<
" events found!" <<
Endl;
694 else if (!bCutMaxDepth)
695 Log() << kVERBOSE <<
"Warning: Maximum depth reached: "
696 << GetMaxDepth() <<
Endl;
698 Log() << kWARNING <<
"<PDEFoam::PeekMax>: no more candidate cells (drivMax>0) found for further splitting." <<
Endl;
719 if(fLastCe+1 >= fNCells)
Log() << kFATAL <<
"Buffer limit is reached, fLastCe=fnBuf" <<
Endl;
726 if( kBest<0 || kBest>=fDim )
Log() << kFATAL <<
"Wrong kBest" <<
Endl;
732 Int_t d1 = CellFill(1, cell);
733 Int_t d2 = CellFill(1, cell);
737 Explore( (fCells[d1]) );
738 Explore( (fCells[d2]) );
753 std::vector<Double_t> xvec;
754 xvec.reserve(GetTotDim());
755 for (
Int_t idim = 0; idim < GetTotDim(); ++idim)
756 xvec.push_back( VarTransformInvers(idim, xRand[idim]) );
758 return GetDistr()->Density(xvec, event_density);
769 fTimer->Init(fNCells);
774 while ( (fLastCe+2) < fNCells ) {
777 if ( (iCell<0) || (iCell>fLastCe) ) {
778 Log() << kVERBOSE <<
"Break: "<< fLastCe+1 <<
" cells created" <<
Endl;
780 for (
Long_t jCell=fLastCe+1; jCell<fNCells; jCell++)
781 delete fCells[jCell];
785 newCell = fCells[iCell];
789 if ( Divide( newCell )==0)
break;
794 Log() << kVERBOSE << GetNActiveCells() <<
" active cells created" <<
Endl;
804 if(fDim==0)
Log() << kFATAL <<
"SetInhiDiv: fDim=0" <<
Endl;
806 fInhiDiv =
new Int_t[ fDim ];
807 for(
Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
810 if( ( 0<=iDim) && (iDim<fDim)) {
811 fInhiDiv[iDim] = inhiDiv;
813 Log() << kFATAL <<
"Wrong iDim" <<
Endl;
824 Int_t errors, warnings;
828 errors = 0; warnings = 0;
829 if (level==1)
Log() << kVERBOSE <<
"Performing consistency checks for created foam" <<
Endl;
830 for(iCell=1; iCell<=fLastCe; iCell++) {
831 cell = fCells[iCell];
836 if (level==1)
Log() << kFATAL <<
"ERROR: Cell's no %d has only one daughter " << iCell <<
Endl;
840 if (level==1)
Log() << kFATAL <<
"ERROR: Cell's no %d has no daughter and is inactive " << iCell <<
Endl;
844 if (level==1)
Log() << kFATAL <<
"ERROR: Cell's no %d has two daughters and is active " << iCell <<
Endl;
848 if( (cell->
GetPare())!=fCells[0] ) {
851 if (level==1)
Log() << kFATAL <<
"ERROR: Cell's no %d parent not pointing to this cell " << iCell <<
Endl;
859 if (level==1)
Log() << kFATAL <<
"ERROR: Cell's no %d daughter 0 not pointing to this cell " << iCell <<
Endl;
865 if (level==1)
Log() << kFATAL <<
"ERROR: Cell's no %d daughter 1 not pointing to this cell " << iCell <<
Endl;
870 if(level==1)
Log() << kFATAL <<
"ERROR: Cell no. " << iCell <<
" has Volume of <1E-50" <<
Endl;
875 for(iCell=0; iCell<=fLastCe; iCell++) {
876 cell = fCells[iCell];
879 if(level==1)
Log() << kFATAL <<
"ERROR: Cell no. " << iCell <<
" is active but Volume is 0 " <<
Endl;
884 Log() << kVERBOSE <<
"Check has found " << errors <<
" errors and " << warnings <<
" warnings." <<
Endl;
887 Info(
"CheckAll",
"Check - found total %d errors \n",errors);
897 if (iCell < 0 || iCell > fLastCe) {
898 Log() << kWARNING <<
"<PrintCell(iCell=" << iCell
899 <<
")>: cell number " << iCell <<
" out of bounds!"
905 fCells[iCell]->GetHcub(cellPosi,cellSize);
906 Int_t kBest = fCells[iCell]->GetBest();
907 Double_t xBest = fCells[iCell]->GetXdiv();
909 Log() <<
"Cell[" << iCell <<
"]={ ";
910 Log() <<
" " << fCells[iCell] <<
" " <<
Endl;
911 Log() <<
" Xdiv[abs. coord.]="
912 << VarTransformInvers(kBest,cellPosi[kBest] + xBest*cellSize[kBest])
914 Log() <<
" Abs. coord. = (";
915 for (
Int_t idim=0; idim<fDim; idim++) {
916 Log() <<
"dim[" << idim <<
"]={"
917 << VarTransformInvers(idim,cellPosi[idim]) <<
","
918 << VarTransformInvers(idim,cellPosi[idim] + cellSize[idim])
924 fCells[iCell]->
Print(
"1");
926 Log() <<
"Elements: [";
930 if (i>0)
Log() <<
", ";
931 Log() << GetCellElement(fCells[iCell], i);
944 for(
Long_t iCell=0; iCell<=fLastCe; iCell++)
958 std::vector<Float_t> values = ev->
GetValues();
959 std::vector<Float_t> tvalues = VarTransform(values);
964 SetCellElement(cell, 0, GetCellElement(cell, 0) + wt);
965 SetCellElement(cell, 1, GetCellElement(cell, 1) + wt*wt);
975 Log() << kVERBOSE <<
"Delete cell elements" <<
Endl;
976 for (
Long_t iCell = 0; iCell < fNCells; ++iCell) {
977 TObject* elements = fCells[iCell]->GetElement();
980 fCells[iCell]->SetElement(NULL);
1018 std::vector<Float_t> txvec(VarTransform(xvec));
1020 return GetCellValue(FindCell(txvec), cv);
1022 return kernel->
Estimate(
this, txvec, cv);
1046 std::map<Int_t,Float_t> txvec;
1047 for (std::map<Int_t,Float_t>::const_iterator it=xvec.begin(); it!=xvec.end(); ++it)
1048 txvec.insert(std::pair<Int_t, Float_t>(it->first, VarTransform(it->first, it->second)));
1051 std::vector<PDEFoamCell*> cells = FindCells(txvec);
1054 std::vector<Float_t> cell_values;
1055 cell_values.reserve(cells.size());
1056 for (std::vector<PDEFoamCell*>::const_iterator cell_it=cells.begin();
1057 cell_it != cells.end(); ++cell_it)
1058 cell_values.push_back(GetCellValue(*cell_it, cv));
1082 PDEFoamVect cellPosi0(GetTotDim()), cellSize0(GetTotDim());
1090 cell0->
GetHcub(cellPosi0,cellSize0);
1092 if (xvec.at(idim)<=cellPosi0[idim]+cellSize0[idim])
1120 PDEFoamVect cellPosi0(GetTotDim()), cellSize0(GetTotDim());
1128 map<Int_t, Float_t>::const_iterator it = txvec.find(idim);
1130 if (it != txvec.end()){
1134 cell0->
GetHcub(cellPosi0,cellSize0);
1136 if (it->second <= cellPosi0[idim] + cellSize0[idim])
1142 FindCells(txvec, cell->
GetDau0(), cells);
1143 FindCells(txvec, cell->
GetDau1(), cells);
1147 cells.push_back(cell);
1168 std::map<Int_t, Float_t> txvec_map;
1169 for (
UInt_t i=0; i<txvec.size(); ++i)
1170 txvec_map.insert(std::pair<Int_t, Float_t>(i, txvec.at(i)));
1173 std::vector<PDEFoamCell*> cells(0);
1176 FindCells(txvec_map, fCells[0], cells);
1199 std::vector<PDEFoamCell*> cells(0);
1202 FindCells(txvec, fCells[0], cells);
1221 if ( GetTotDim()!=1 )
1222 Log() << kFATAL <<
"<Draw1Dim>: function can only be used for 1-dimensional foams!"
1228 h1=
new TH1D(hname,
"1-dimensional Foam", nbin, fXmin[0], fXmax[0]);
1230 if (!
h1)
Log() << kFATAL <<
"ERROR: Can not create histo" << hname <<
Endl;
1235 std::vector<Float_t> txvec;
1238 if (kernel != NULL) {
1240 val = kernel->
Estimate(
this, txvec, cell_value);
1242 val = GetCellValue(FindCell(txvec), cell_value);
1273 if ((idim1>=GetTotDim()) || (idim1<0) ||
1274 (idim2>=GetTotDim()) || (idim2<0) ||
1276 Log() << kFATAL <<
"<Project2>: wrong dimensions given: "
1277 << idim1 <<
", " << idim2 <<
Endl;
1283 Log() << kWARNING <<
"Warning: number of bins too big: " << nbin
1284 <<
" Using 1000 bins for each dimension instead." <<
Endl;
1286 }
else if (nbin<1) {
1287 Log() << kWARNING <<
"Wrong bin number: " << nbin
1288 <<
"; set nbin=50" <<
Endl;
1298 h1=
new TH2D(hname.
Data(),
Form(
"var%d vs var%d",idim1,idim2), nbin, fXmin[idim1], fXmax[idim1], nbin, fXmin[idim2], fXmax[idim2]);
1300 if (!
h1)
Log() << kFATAL <<
"ERROR: Can not create histo" << hname <<
Endl;
1308 std::map<Int_t, Float_t> txvec;
1314 std::vector<TMVA::PDEFoamCell*> cells = FindCells(txvec);
1319 for (std::vector<TMVA::PDEFoamCell*>::const_iterator it = cells.begin();
1320 it != cells.end(); ++it) {
1322 PDEFoamVect cellPosi(GetTotDim()), cellSize(GetTotDim());
1323 (*it)->GetHcub(cellPosi,cellSize);
1326 std::vector<Float_t> tvec;
1327 for (
Int_t i=0; i<GetTotDim(); ++i) {
1328 if ( i != idim1 && i != idim2 )
1329 tvec.push_back(cellPosi[i] + 0.5*cellSize[i]);
1331 tvec.push_back(txvec[i]);
1333 if (kernel != NULL) {
1335 sum_cv += kernel->
Estimate(
this, tvec, cell_value);
1337 sum_cv += GetCellValue(FindCell(tvec), cell_value);
1362 return GetCellElement(cell, 0);
1365 return GetCellElement(cell, 1);
1367 case kValueDensity: {
1371 return GetCellValue(cell,
kValue)/volume;
1375 Log() << kWARNING <<
"<GetCellDensity(cell)>: ERROR: cell volume"
1376 <<
" negative or zero!"
1377 <<
" ==> return cell density 0!"
1378 <<
" cell volume=" << volume
1379 <<
" cell entries=" << GetCellValue(cell,
kValue) <<
Endl;
1381 Log() << kWARNING <<
"<GetCellDensity(cell)>: WARNING: cell volume"
1382 <<
" close to zero!"
1383 <<
" cell volume: " << volume <<
Endl;
1405 Log() << kFATAL <<
"<GetCellValue>: unknown cell value" <<
Endl;
1447 Log() << kFATAL <<
"<SetCellElement> ERROR: cell element is not a TVectorD*" <<
Endl;
1463 Log() << kINFO <<
"Elapsed time: " + fTimer->GetElapsedTime()
1470 if (fNCells >= 100) modulo =
Int_t(fNCells/100);
1471 if (fLastCe%modulo == 0) fTimer->DrawProgressBar( fLastCe );
1514 if (GetTotDim() != 2)
1515 Log() << kFATAL <<
"RootPlot2dim() can only be used with "
1516 <<
"two-dimensional foams!" <<
Endl;
1519 ECellValue cell_value =
kValue;
1524 }
else if (opt.
Contains(
"rms_ov_mean")){
1525 cell_value = kRmsOvMean;
1532 plotcellnumber =
kTRUE;
1535 std::ofstream outfile(filename, std::ios::out);
1537 outfile<<
"{" << std::endl;
1542 outfile <<
"TColor *graycolors[100];" << std::endl;
1543 outfile <<
"for (Int_t i=0.; i<100; i++)" << std::endl;
1544 outfile <<
" graycolors[i]=new TColor(1000+i, 1-(Float_t)i/100.,1-(Float_t)i/100.,1-(Float_t)i/100.);"<< std::endl;
1547 outfile <<
"cMap = new TCanvas(\"" << fName <<
"\",\"Cell Map for "
1548 << fName <<
"\",600,600);" << std::endl;
1550 outfile<<
"TBox*a=new TBox();"<<std::endl;
1551 outfile<<
"a->SetFillStyle(0);"<<std::endl;
1552 outfile<<
"a->SetLineWidth(4);"<<std::endl;
1553 outfile<<
"TBox *b1=new TBox();"<<std::endl;
1554 outfile<<
"TText*t=new TText();"<<std::endl;
1556 outfile << (
colors ?
"gStyle->SetPalette(1, 0);" :
"gStyle->SetPalette(0);")
1558 outfile <<
"b1->SetFillStyle(1001);"<<std::endl;
1559 outfile<<
"TBox *b2=new TBox();"<<std::endl;
1560 outfile <<
"b2->SetFillStyle(0);"<<std::endl;
1563 outfile <<
"b1->SetFillStyle(0);"<<std::endl;
1575 for (
Long_t iCell=1; iCell<=fLastCe; iCell++) {
1576 if ( fCells[iCell]->GetStat() == 1) {
1577 Float_t value = GetCellValue(fCells[iCell], cell_value);
1584 outfile <<
"// observed minimum and maximum of distribution: " << std::endl;
1585 outfile <<
"// Float_t zmin = "<< zmin <<
";" << std::endl;
1586 outfile <<
"// Float_t zmax = "<< zmax <<
";" << std::endl;
1589 outfile <<
"// used minimum and maximum of distribution (taking into account log scale if applicable): " << std::endl;
1590 outfile <<
"Float_t zmin = "<< zmin <<
";" << std::endl;
1591 outfile <<
"Float_t zmax = "<< zmax <<
";" << std::endl;
1597 Float_t scale = (ncolors-1)/(zmax - zmin);
1598 PDEFoamVect cellPosi(GetTotDim()), cellSize(GetTotDim());
1602 outfile <<
"// =========== Rectangular cells ==========="<< std::endl;
1603 for (
Long_t iCell=1; iCell<=fLastCe; iCell++) {
1604 if ( fCells[iCell]->GetStat() == 1) {
1605 fCells[iCell]->GetHcub(cellPosi,cellSize);
1606 x1 = offs+lpag*(cellPosi[0]);
1607 y1 = offs+lpag*(cellPosi[1]);
1608 x2 = offs+lpag*(cellPosi[0]+cellSize[0]);
1609 y2 = offs+lpag*(cellPosi[1]+cellSize[1]);
1613 Float_t value = GetCellValue(fCells[iCell], cell_value);
1620 color = 1000+(
Int_t((value-zmin)*scale));
1623 outfile <<
"b1->SetFillColor(" << color <<
");" << std::endl;
1627 outfile<<
"b1->DrawBox("<<
x1<<
","<<y1<<
","<<
x2<<
","<<y2<<
");"<<std::endl;
1629 outfile<<
"b2->DrawBox("<<
x1<<
","<<y1<<
","<<
x2<<
","<<y2<<
");"<<std::endl;
1632 if (plotcellnumber) {
1633 outfile<<
"t->SetTextColor(4);"<<std::endl;
1635 outfile<<
"t->SetTextSize(0.025);"<<std::endl;
1636 else if(fLastCe<251)
1637 outfile<<
"t->SetTextSize(0.015);"<<std::endl;
1639 outfile<<
"t->SetTextSize(0.008);"<<std::endl;
1640 x = offs+lpag*(cellPosi[0]+0.5*cellSize[0]);
1641 y = offs+lpag*(cellPosi[1]+0.5*cellSize[1]);
1642 outfile<<
"t->DrawText("<<
x<<
","<<
y<<
","<<
"\""<<iCell<<
"\""<<
");"<<std::endl;
1646 outfile<<
"// ============== End Rectangles ==========="<< std::endl;
1648 outfile <<
"}" << std::endl;
1659 GetDistr()->FillBinarySearchTree(ev);
1668 if(fDistr)
delete fDistr;