STAT/TKDNodeInfo.cxx

   1 #include "TKDNodeInfo.h"
   2
   3 #include "TVectorD.h"
   4 #include "TMatrixD.h"
   5 #include "TMath.h"
   6
   7 ClassImp(TKDNodeInfo)
   8
   9
  10 //_________________________________________________________________
  11 TKDNodeInfo::TKDNodeInfo(const Int_t dim):
  12         TObject()
  13         ,fNDim(3*dim)
  14         ,fData(0x0)
  15         ,fCov(0x0)
  16         ,fPar(0x0)
  17 {
  18         fVal[0] = 0.; fVal[1] = 0.;
  19         Build(dim);
  20 }
  21
  22 //_________________________________________________________________
  23 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
  24         TObject((TObject&) ref)
  25         ,fNDim(fNDim)
  26         ,fData(0x0)
  27         ,fCov(0x0)
  28         ,fPar(0x0)
  29 {
  30         Build(fNDim/3);
  31
  32         memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
  33         fVal[0] = ref.fVal[0];
  34         fVal[1] = ref.fVal[1];
  35         if(ref.fCov) (*fCov) = (*ref.fCov);
  36         if(ref.fPar) (*fPar) = (*ref.fPar);
  37 }
  38
  39
  40 //_________________________________________________________________
  41 TKDNodeInfo::~TKDNodeInfo()
  42 {
  43         if(fData) delete [] fData;
  44         if(fCov){
  45                 delete fPar;
  46                 delete fCov;
  47         }
  48 }
  49
  50 //_________________________________________________________________
  51 TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
  52 {
  53 //      Info("operator==()", "...");
  54
  55         Int_t ndim = fNDim/3;
  56         if(fNDim != ref.fNDim){
  57                 fNDim = ref.fNDim;
  58                 Build(ndim);
  59         }
  60         memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
  61         fVal[0] = ref.fVal[0];
  62         fVal[1] = ref.fVal[1];
  63         if(ref.fCov) (*fCov) = (*ref.fCov);
  64         if(ref.fPar) (*fPar) = (*ref.fPar);
  65
  66         return *this;
  67 }
  68
  69 //_________________________________________________________________
  70 void TKDNodeInfo::Build(const Int_t dim)
  71 {
  72 // Allocate/Reallocate space for this node.
  73
  74 //      Info("Build()", "...");
  75
  76         if(!dim) return;
  77
  78         Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
  79         if(fData) delete [] fData;
  80         fData = new Float_t[fNDim];
  81         if(fCov){
  82                 fCov->ResizeTo(lambda, lambda);
  83                 fPar->ResizeTo(lambda);
  84         }
  85         return;
  86 }
  87
  88 //_________________________________________________________________
  89 void TKDNodeInfo::Print()
  90 {
  91         Int_t dim = fNDim/3;
  92         printf("x[");
  93         for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
  94         printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
  95
  96         Float_t *bounds = &fData[dim];
  97         printf("range[");
  98         for(int idim=0; idim<dim; idim++) printf("(%f %f) ", fData[2*idim], fData[2*idim+1]);
  99         printf("]\n");
 100
 101         printf("Fit parameters : ");
 102         if(!fCov){
 103                 printf("Not defined.\n");
 104                 return;
 105         }
 106
 107         Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
 108         for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
 109         printf("\n");
 110 }
 111
 112 //_________________________________________________________________
 113 void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
 114 {
 115         if(!fCov){
 116                 fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
 117                 fPar = new TVectorD(par.GetNrows());
 118         }
 119         (*fPar) = par;
 120         (*fCov) = cov;
 121 }
 122
 123 //_________________________________________________________________
 124 Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
 125 {
 126 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
 127
 128         Int_t ndim = fNDim/3;
 129         if(ndim>10) return 0.; // support only up to 10 dimensions
 130
 131         Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
 132         Double_t fdfdp[66];
 133         Int_t ipar = 0;
 134         fdfdp[ipar++] = 1.;
 135         for(int idim=0; idim<ndim; idim++){
 136                 fdfdp[ipar++] = point[idim];
 137                 for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
 138         }
 139
 140         // calculate estimation
 141         result =0.; error = 0.;
 142         for(int i=0; i<lambda; i++){
 143                 result += fdfdp[i]*(*fPar)(i);
 144                 for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
 145         }
 146         error = TMath::Sqrt(error);
 147
 148         //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
 149
 150         return 0.;
 151 }
 152