STAT/TKDNodeInfo.cxx

   1 ////////////////////////////////////////////////////////
   2 //
   3 // Bucket representation for TKDInterpolator(Base)
   4 //
   5 // The class store data and provides the interface to draw and print.
   6 // The bucket - generalized histogram bin in N dimensions is represented by unprocessed data like
   7 //   - experimental PDF value and statistical error
   8 //   - COG position (n-tuplu)
   9 //   - boundaries
  10 // and interpolated data like
  11 //   - parameters of the local parabolic fit
  12 //   - their covariance matrix
  13 //
  14 // For drawing 2D projections the helper class TKDNodeInfo::TKDNodeDraw is used.
  15 //
  16 // Author Alexandru Bercuci <A.Bercuci@gsi.de>
  17 //
  18 ////////////////////////////////////////////////////////
  19
  20 #include "TKDNodeInfo.h"
  21
  22 #include "TRandom.h"
  23 #include "TVectorD.h"
  24 #include "TMatrixD.h"
  25 #include "TMath.h"
  26
  27 ClassImp(TKDNodeInfo)
  28 ClassImp(TKDNodeInfo::TKDNodeDraw)
  29
  30
  31 //_________________________________________________________________
  32 TKDNodeInfo::TKDNodeInfo(Int_t dim):
  33   TObject()
  34   ,fNDim(3*dim)
  35   ,fData(0x0)
  36   ,fCov(0x0)
  37   ,fPar(0x0)
  38 {
  39   // Default constructor
  40   fVal[0] = 0.; fVal[1] = 0.;
  41   Build(dim);
  42 }
  43
  44 //_________________________________________________________________
  45 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
  46   TObject((TObject&) ref)
  47   ,fNDim(fNDim)
  48   ,fData(0x0)
  49   ,fCov(0x0)
  50   ,fPar(0x0)
  51 {
  52   // Copy constructor
  53   Build(fNDim/3);
  54
  55   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
  56   fVal[0] = ref.fVal[0];
  57   fVal[1] = ref.fVal[1];
  58   if(ref.fCov) (*fCov) = (*ref.fCov);
  59   if(ref.fPar) (*fPar) = (*ref.fPar);
  60 }
  61
  62
  63 //_________________________________________________________________
  64 TKDNodeInfo::~TKDNodeInfo()
  65 {
  66   // Destructor
  67   if(fData) delete [] fData;
  68   if(fCov){
  69     delete fPar;
  70     delete fCov;
  71   }
  72 }
  73
  74 //_________________________________________________________________
  75 TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
  76 {
  77 //      Info("operator==()", "...");
  78
  79   Int_t ndim = fNDim/3;
  80   if(fNDim != ref.fNDim){
  81     fNDim = ref.fNDim;
  82     Build(ndim);
  83   }
  84   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
  85   fVal[0] = ref.fVal[0];
  86   fVal[1] = ref.fVal[1];
  87   if(ref.fCov) (*fCov) = (*ref.fCov);
  88   if(ref.fPar) (*fPar) = (*ref.fPar);
  89
  90   return *this;
  91 }
  92
  93 //_________________________________________________________________
  94 void TKDNodeInfo::Build(Int_t dim)
  95 {
  96 // Allocate/Reallocate space for this node.
  97
  98 //      Info("Build()", "...");
  99
 100   if(!dim) return;
 101   fNDim = 3*dim;
 102   Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
 103   if(fData) delete [] fData;
 104   fData = new Float_t[fNDim];
 105   if(fCov){
 106     fCov->ResizeTo(lambda, lambda);
 107     fPar->ResizeTo(lambda);
 108   }
 109   return;
 110 }
 111
 112 //_________________________________________________________________
 113 void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
 114 {
 115   Build(ndim);
 116   memcpy(fData, data, fNDim*sizeof(Float_t));
 117   fVal[0]=pdf[0]; fVal[1]=pdf[1];
 118 }
 119
 120
 121 //_________________________________________________________________
 122 void TKDNodeInfo::Print(const Option_t *) const
 123 {
 124   // Print the content of the node
 125   Int_t dim = fNDim/3;
 126   printf("x[");
 127   for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
 128   printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
 129
 130   Float_t *bounds = &fData[dim];
 131   printf("range[");
 132   for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
 133   printf("]\n");
 134
 135   printf("Fit parameters : ");
 136   if(!fCov){
 137     printf("Not defined.\n");
 138     return;
 139   }
 140
 141   //    Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
 142   for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
 143   printf("\n");
 144 }
 145
 146 //_________________________________________________________________
 147 void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
 148 {
 149   // Store the parameters and the covariance in the node
 150   if(!fCov){
 151     fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
 152     fPar = new TVectorD(par.GetNrows());
 153   }
 154   (*fPar) = par;
 155   (*fCov) = cov;
 156 }
 157
 158 //_________________________________________________________________
 159 Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
 160 {
 161 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
 162
 163   Int_t ndim = fNDim/3;
 164   if(ndim>10) return 0.; // support only up to 10 dimensions
 165
 166   Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
 167   Double_t fdfdp[66];
 168   Int_t ipar = 0;
 169   fdfdp[ipar++] = 1.;
 170   for(int idim=0; idim<ndim; idim++){
 171     fdfdp[ipar++] = point[idim];
 172     for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
 173   }
 174
 175   // calculate estimation
 176   result =0.; error = 0.;
 177   for(int i=0; i<lambda; i++){
 178     result += fdfdp[i]*(*fPar)(i);
 179     for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
 180   }
 181   error = TMath::Sqrt(error);
 182
 183   //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
 184
 185   return 0.;
 186 }
 187
 188
 189
 190 //_________________________________________________________________
 191 TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
 192   :TBox()
 193   ,fCOG()
 194   ,fNode(0x0)
 195 {
 196   SetFillStyle(3002);
 197   SetFillColor(50+Int_t(gRandom->Uniform()*50.));
 198
 199   fCOG.SetMarkerStyle(3);
 200   fCOG.SetMarkerSize(.7);
 201   fCOG.SetMarkerColor(2);
 202 }
 203
 204
 205 //_________________________________________________________________
 206 void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
 207 {
 208   TBox::Draw(option);
 209   fCOG.Draw("p");
 210 }
 211
 212 //_________________________________________________________________
 213 void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
 214 {
 215   fNode=node;
 216   const Float_t kBorder = 0.;//1.E-4;
 217   Float_t *bounds = &(node->Data()[size]);
 218   fX1=bounds[2*ax1]+kBorder;
 219   fX2=bounds[2*ax1+1]-kBorder;
 220   fY1=bounds[2*ax2]+kBorder;
 221   fY2=bounds[2*ax2+1]-kBorder;
 222
 223   Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
 224   fCOG.SetX(x); fCOG.SetY(y);
 225 }
 226
 227
 228 //_________________________________________________________________
 229 void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
 230 {
 231   if(!fNode) return;
 232   fNode->Print(option);
 233 }