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