+////////////////////////////////////////////////////////
+//
+// Bucket representation for TKDInterpolator(Base)
+//
+// The class store data and provides the interface to draw and print.
+// The bucket - generalized histogram bin in N dimensions is represented by unprocessed data like
+// - experimental PDF value and statistical error
+// - COG position (n-tuplu)
+// - boundaries
+// and interpolated data like
+// - parameters of the local parabolic fit
+// - their covariance matrix
+//
+// For drawing 2D projections the helper class TKDNodeInfo::TKDNodeDraw is used.
+//
+// Author Alexandru Bercuci <A.Bercuci@gsi.de>
+//
+////////////////////////////////////////////////////////
+
#include "TKDNodeInfo.h"
-#include "TRandom.h"
#include "TVectorD.h"
#include "TMatrixD.h"
+#include "TRandom.h"
#include "TMath.h"
ClassImp(TKDNodeInfo)
TKDNodeInfo::TKDNodeInfo(Int_t dim):
TObject()
,fNDim(3*dim)
- ,fData(0x0)
- ,fCov(0x0)
- ,fPar(0x0)
+ ,fData(NULL)
+ ,fNpar(0)
+ ,fNcov(0)
+ ,fPar(NULL)
+ ,fCov(NULL)
{
// Default constructor
fVal[0] = 0.; fVal[1] = 0.;
//_________________________________________________________________
TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
TObject((TObject&) ref)
- ,fNDim(fNDim)
- ,fData(0x0)
- ,fCov(0x0)
- ,fPar(0x0)
+ ,fNDim(ref.fNDim)
+ ,fData(NULL)
+ ,fNpar(0)
+ ,fNcov(0)
+ ,fPar(NULL)
+ ,fCov(NULL)
{
// Copy constructor
Build(fNDim/3);
+ fData = new Float_t[fNDim];
memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
fVal[0] = ref.fVal[0];
fVal[1] = ref.fVal[1];
- if(ref.fCov) (*fCov) = (*ref.fCov);
- if(ref.fPar) (*fPar) = (*ref.fPar);
+ if(ref.fNpar&&ref.fPar){
+ fNpar = ref.fNpar;
+ fPar=new Double_t[fNpar];
+ memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
+ }
+ if(ref.fNcov && ref.fCov){
+ fNcov = ref.fNcov;
+ fCov=new Double_t[fNcov];
+ memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
+ }
}
{
// Destructor
if(fData) delete [] fData;
- if(fCov){
- delete fPar;
- delete fCov;
- }
+ if(fPar) delete [] fPar;
+ if(fCov) delete [] fCov;
}
//_________________________________________________________________
{
// Info("operator==()", "...");
+ if(this == &ref) return *this;
Int_t ndim = fNDim/3;
if(fNDim != ref.fNDim){
fNDim = ref.fNDim;
memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
fVal[0] = ref.fVal[0];
fVal[1] = ref.fVal[1];
- if(ref.fCov) (*fCov) = (*ref.fCov);
- if(ref.fPar) (*fPar) = (*ref.fPar);
-
+ if(ref.fNpar&&ref.fPar){
+ fNpar = ref.fNpar;
+ fPar=new Double_t[fNpar];
+ memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
+ }
+ if(ref.fNcov && ref.fCov){
+ fNcov = ref.fNcov;
+ fCov=new Double_t[fNcov];
+ memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
+ }
return *this;
}
// Info("Build()", "...");
if(!dim) return;
-
- Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
+ fNDim = 3*dim;
if(fData) delete [] fData;
fData = new Float_t[fNDim];
- if(fCov){
- fCov->ResizeTo(lambda, lambda);
- fPar->ResizeTo(lambda);
- }
return;
}
//_________________________________________________________________
-void TKDNodeInfo::Print(const Option_t *) const
+void TKDNodeInfo::Bootstrap()
+{
+ if(!fNpar || !fPar) return;
+
+ Int_t ndim = Int_t(.5*(TMath::Sqrt(1.+8.*fNpar)-1.))-1;
+ fNDim = 3*ndim;
+}
+
+//_________________________________________________________________
+void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
+{
+ Build(ndim);
+ memcpy(fData, data, fNDim*sizeof(Float_t));
+ fVal[0]=pdf[0]; fVal[1]=pdf[1];
+}
+
+
+//_________________________________________________________________
+void TKDNodeInfo::Print(const Option_t *opt) const
{
// Print the content of the node
- Int_t dim = fNDim/3;
+ Int_t dim = Int_t(fNDim/3.);
printf("x[");
- for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
+ for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
- // Float_t *bounds = &fData[dim];
- printf("range[");
- for(int idim=0; idim<dim; idim++) printf("(%f %f) ", fData[2*idim], fData[2*idim+1]);
- printf("]\n");
-
- printf("Fit parameters : ");
- if(!fCov){
- printf("Not defined.\n");
- return;
+ if(fData){
+ Float_t *bounds = &fData[dim];
+ printf("range[");
+ for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
+ printf("]\n");
+ }
+ if(strcmp(opt, "a")!=0) return;
+
+ if(fNpar){
+ printf("Fit parameters : \n");
+ for(int ip=0; ip<fNpar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
+ printf("\n");
+ }
+ if(!fNcov) return;
+ for(int ip(0), n(0); ip<fNpar; ip++){
+ for(int jp(ip); jp<fNpar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
+ printf("\n");
}
-
- // Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
- for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
- printf("\n");
}
//_________________________________________________________________
-void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
+void TKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
{
- // Store the parameters and the covariance in the node
- if(!fCov){
- fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
- fPar = new TVectorD(par.GetNrows());
- }
- (*fPar) = par;
- (*fCov) = cov;
+// Store the parameters and the covariance in the node
+
+ if(!fPar){SetNpar(); fPar = new Double_t[fNpar];}
+ for(int ip=0; ip<fNpar; ip++) fPar[ip] = (*par)[ip];
+
+ if(!cov) return;
+ if(!fCov){SetNcov(); fCov = new Double_t[fNcov];}
+ for(int ip(0), np(0); ip<fNpar; ip++)
+ for(int jp=ip; jp<fNpar; jp++) fCov[np++] = (*cov)(ip,jp);
}
//_________________________________________________________________
-Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
+Bool_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
{
// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
- Int_t ndim = fNDim/3;
- if(ndim>10) return 0.; // support only up to 10 dimensions
+ Int_t ndim = Int_t(fNDim/3.);
+ if(ndim>10) return kFALSE; // support only up to 10 dimensions
+ //printf("ndim[%d] npar[%d] ncov[%d]\n", ndim, fNpar, fNcov);
- Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
- Double_t fdfdp[66];
+ Double_t fdfdp[66]; memset(fdfdp, 0, ndim*sizeof(Double_t));
Int_t ipar = 0;
fdfdp[ipar++] = 1.;
for(int idim=0; idim<ndim; idim++){
// calculate estimation
result =0.; error = 0.;
- for(int i=0; i<lambda; i++){
- result += fdfdp[i]*(*fPar)(i);
- for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
+ for(int i=0; i<fNpar; i++) result += fdfdp[i]*fPar[i];
+ if(!fNcov) return kTRUE;
+
+ for(int i(0), n(0); i<fNpar; i++){
+ error += fdfdp[i]*fdfdp[i]*fCov[n++];
+ for(int j(i+1); j<fNpar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
}
error = TMath::Sqrt(error);
//printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
- return 0.;
+ return kTRUE;
}
//_________________________________________________________________
-TKDNodeInfo::TKDNodeDraw::TKDNodeDraw() : TBox(), fCOG()
+TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
+ :TBox()
+ ,fCOG()
+ ,fNode(NULL)
{
SetFillStyle(3002);
SetFillColor(50+Int_t(gRandom->Uniform()*50.));
//_________________________________________________________________
void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
{
+ fNode=node;
const Float_t kBorder = 0.;//1.E-4;
Float_t *bounds = &(node->Data()[size]);
- //printf("x1[%f] x2[%f] y1[%f] y2[%f]\n", bounds[2*ax1], bounds[2*ax1+1], bounds[2*ax2], bounds[2*ax2+1]);
fX1=bounds[2*ax1]+kBorder;
fX2=bounds[2*ax1+1]-kBorder;
fY1=bounds[2*ax2]+kBorder;
fY2=bounds[2*ax2+1]-kBorder;
Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
- //printf("x[%f] y[%f]\n", x, y);
fCOG.SetX(x); fCOG.SetY(y);
}
+//_________________________________________________________________
+void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
+{
+ if(!fNode) return;
+ fNode->Print(option);
+}