]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STAT/TKDNodeInfo.cxx
From Michael Knichel: The THnSparses are not merged anymore, instead a set of 1D...
[u/mrichter/AliRoot.git] / STAT / TKDNodeInfo.cxx
index 85b25ebe6d7da57ea1bc1aa436b82ec3b879a9b6..915a688183b1d159cd91ed49ba14688d307f6d4f 100644 (file)
@@ -19,9 +19,9 @@
 
 #include "TKDNodeInfo.h"
 
-#include "TRandom.h"
 #include "TVectorD.h"
 #include "TMatrixD.h"
+#include "TRandom.h"
 #include "TMath.h"
 
 ClassImp(TKDNodeInfo)
@@ -32,9 +32,11 @@ ClassImp(TKDNodeInfo::TKDNodeDraw)
 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.;
@@ -45,9 +47,11 @@ TKDNodeInfo::TKDNodeInfo(Int_t dim):
 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
   TObject((TObject&) ref)
   ,fNDim(fNDim)
-  ,fData(0x0)
-  ,fCov(0x0)
-  ,fPar(0x0)
+  ,fData(NULL)
+  ,fNpar(0)
+  ,fNcov(0)
+  ,fPar(NULL)
+  ,fCov(NULL)
 {
   // Copy constructor
   Build(fNDim/3);
@@ -55,8 +59,16 @@ TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
   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));
+  }
 }
 
 
@@ -65,10 +77,8 @@ TKDNodeInfo::~TKDNodeInfo()
 {
   // Destructor
   if(fData) delete [] fData;
-  if(fCov){
-    delete fPar;
-    delete fCov;
-  }
+  if(fPar) delete [] fPar;
+  if(fCov) delete [] fCov;
 }
 
 //_________________________________________________________________
@@ -84,9 +94,16 @@ TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
   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;
 }
 
@@ -98,17 +115,21 @@ void TKDNodeInfo::Build(Int_t dim)
 //     Info("Build()", "...");
 
   if(!dim) return;
-  fNDim = 3*dim;  
-  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::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)
 {
@@ -119,52 +140,58 @@ void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
 
 
 //_________________________________________________________________
-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;
+  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) ", bounds[2*idim], bounds[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++){
@@ -174,15 +201,18 @@ Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t
 
   // 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;
 }
 
 
@@ -191,7 +221,7 @@ Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t
 TKDNodeInfo::TKDNodeDraw::TKDNodeDraw() 
   :TBox()
   ,fCOG()
-  ,fNode(0x0)
+  ,fNode(NULL)
 {
   SetFillStyle(3002);
   SetFillColor(50+Int_t(gRandom->Uniform()*50.));