+////////////////////////////////////////////////////////
+//
+// 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 "TVectorD.h"
#include "TMatrixD.h"
+#include "TRandom.h"
#include "TMath.h"
ClassImp(TKDNodeInfo)
+ClassImp(TKDNodeInfo::TKDNodeDraw)
//_________________________________________________________________
TKDNodeInfo::TKDNodeInfo(Int_t dim):
- TObject()
- ,fNDim(3*dim)
- ,fData(0x0)
- ,fCov(0x0)
- ,fPar(0x0)
+ TObject()
+ ,fNDim(3*dim)
+ ,fData(NULL)
+ ,fNpar(0)
+ ,fNcov(0)
+ ,fPar(NULL)
+ ,fCov(NULL)
{
// Default constructor
- fVal[0] = 0.; fVal[1] = 0.;
- Build(dim);
+ fVal[0] = 0.; fVal[1] = 0.;
+ Build(dim);
}
//_________________________________________________________________
TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
- TObject((TObject&) ref)
- ,fNDim(fNDim)
- ,fData(0x0)
- ,fCov(0x0)
- ,fPar(0x0)
+ TObject((TObject&) ref)
+ ,fNDim(ref.fNDim)
+ ,fData(NULL)
+ ,fNpar(0)
+ ,fNcov(0)
+ ,fPar(NULL)
+ ,fCov(NULL)
{
// Copy constructor
- Build(fNDim/3);
+ Build(fNDim/3);
- 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);
+ 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.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));
+ }
}
TKDNodeInfo::~TKDNodeInfo()
{
// Destructor
- if(fData) delete [] fData;
- if(fCov){
- delete fPar;
- delete fCov;
- }
+ if(fData) delete [] fData;
+ if(fPar) delete [] fPar;
+ if(fCov) delete [] fCov;
}
//_________________________________________________________________
TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
{
// Info("operator==()", "...");
-
- Int_t ndim = fNDim/3;
- if(fNDim != ref.fNDim){
- fNDim = ref.fNDim;
- Build(ndim);
- }
- 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);
-
- return *this;
+
+ if(this == &ref) return *this;
+ Int_t ndim = fNDim/3;
+ if(fNDim != ref.fNDim){
+ fNDim = ref.fNDim;
+ Build(ndim);
+ }
+ memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+ fVal[0] = ref.fVal[0];
+ fVal[1] = ref.fVal[1];
+ 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));
- if(fData) delete [] fData;
- fData = new Float_t[fNDim];
- if(fCov){
- fCov->ResizeTo(lambda, lambda);
- fPar->ResizeTo(lambda);
- }
- return;
+ if(!dim) return;
+ fNDim = 3*dim;
+ if(fData) delete [] fData;
+ fData = new Float_t[fNDim];
+ return;
+}
+
+//_________________________________________________________________
+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 *) const
+void TKDNodeInfo::Print(const Option_t *opt) const
{
// Print the content of the node
- Int_t dim = fNDim/3;
- printf("x[");
- for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
- 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;
- }
-
- // 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");
+ Int_t dim = Int_t(fNDim/3.);
+ printf("x[");
+ for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
+ printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
+
+ 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");
+ }
}
//_________________________________________________________________
-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 lambda = 1 + ndim + (ndim*(ndim+1)>>1);
- Double_t fdfdp[66];
- Int_t ipar = 0;
- fdfdp[ipar++] = 1.;
- for(int idim=0; idim<ndim; idim++){
- fdfdp[ipar++] = point[idim];
- for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
- }
-
- // 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);
- }
- error = TMath::Sqrt(error);
-
- //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
-
- return 0.;
+ 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);
+
+ 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++){
+ fdfdp[ipar++] = point[idim];
+ for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
+ }
+
+ // calculate estimation
+ result =0.; error = 0.;
+ 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 kTRUE;
}
+
+
+//_________________________________________________________________
+TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
+ :TBox()
+ ,fCOG()
+ ,fNode(NULL)
+{
+ SetFillStyle(3002);
+ SetFillColor(50+Int_t(gRandom->Uniform()*50.));
+
+ fCOG.SetMarkerStyle(3);
+ fCOG.SetMarkerSize(.7);
+ fCOG.SetMarkerColor(2);
+}
+
+
+//_________________________________________________________________
+void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
+{
+ TBox::Draw(option);
+ fCOG.Draw("p");
+}
+
+//_________________________________________________________________
+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]);
+ 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]);
+ fCOG.SetX(x); fCOG.SetY(y);
+}
+
+
+//_________________________________________________________________
+void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
+{
+ if(!fNode) return;
+ fNode->Print(option);
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff