#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(NULL)
,fNpar(0)
,fNcov(0)
,fPar(NULL)
,fCov(NULL)
{
fVal[0] = 0.; fVal[1] = 0.;
Build(dim);
}
TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
TObject((TObject&) ref)
,fNDim(ref.fNDim)
,fData(NULL)
,fNpar(0)
,fNcov(0)
,fPar(NULL)
,fCov(NULL)
{
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.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()
{
if(fData) delete [] fData;
if(fPar) delete [] fPar;
if(fCov) delete [] fCov;
}
TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
{
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;
}
void TKDNodeInfo::Build(Int_t dim)
{
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 *opt) const
{
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(TVectorD const *par, TMatrixD const *cov)
{
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);
}
Bool_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
{
Int_t ndim = Int_t(fNDim/3.);
if(ndim>10) return kFALSE;
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];
}
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);
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.;
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);
}