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 "TVectorD.h"
  23 #include "TMatrixD.h"
  24 #include "TRandom.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(NULL)
  36   ,fNpar(0)
  37   ,fNcov(0)
  38   ,fPar(NULL)
  39   ,fCov(NULL)
  40 {
  41   // Default constructor
  42   fVal[0] = 0.; fVal[1] = 0.;
  43   Build(dim);
  44 }
  45
  46 //_________________________________________________________________
  47 TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
  48   TObject((TObject&) ref)
  49   ,fNDim(ref.fNDim)
  50   ,fData(NULL)
  51   ,fNpar(0)
  52   ,fNcov(0)
  53   ,fPar(NULL)
  54   ,fCov(NULL)
  55 {
  56   // Copy constructor
  57   Build(fNDim/3);
  58
  59   fData = new Float_t[fNDim];
  60   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
  61   fVal[0] = ref.fVal[0];
  62   fVal[1] = ref.fVal[1];
  63   if(ref.fNpar&&ref.fPar){
  64     fNpar = ref.fNpar;
  65     fPar=new Double_t[fNpar];
  66     memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
  67   }
  68   if(ref.fNcov && ref.fCov){
  69     fNcov = ref.fNcov;
  70     fCov=new Double_t[fNcov];
  71     memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
  72   }
  73 }
  74
  75
  76 //_________________________________________________________________
  77 TKDNodeInfo::~TKDNodeInfo()
  78 {
  79   // Destructor
  80   if(fData) delete [] fData;
  81   if(fPar) delete [] fPar;
  82   if(fCov) delete [] fCov;
  83 }
  84
  85 //_________________________________________________________________
  86 TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
  87 {
  88 //      Info("operator==()", "...");
  89
  90   if(this == &ref) return *this;
  91   Int_t ndim = fNDim/3;
  92   if(fNDim != ref.fNDim){
  93     fNDim = ref.fNDim;
  94     Build(ndim);
  95   }
  96   memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
  97   fVal[0] = ref.fVal[0];
  98   fVal[1] = ref.fVal[1];
  99   if(ref.fNpar&&ref.fPar){
 100     fNpar = ref.fNpar;
 101     fPar=new Double_t[fNpar];
 102     memcpy(fPar, ref.fPar, fNpar*sizeof(Double_t));
 103   }
 104   if(ref.fNcov && ref.fCov){
 105     fNcov = ref.fNcov;
 106     fCov=new Double_t[fNcov];
 107     memcpy(fCov, ref.fCov, fNcov*sizeof(Double_t));
 108   }
 109   return *this;
 110 }
 111
 112 //_________________________________________________________________
 113 void TKDNodeInfo::Build(Int_t dim)
 114 {
 115 // Allocate/Reallocate space for this node.
 116
 117 //      Info("Build()", "...");
 118
 119   if(!dim) return;
 120   fNDim = 3*dim;
 121   if(fData) delete [] fData;
 122   fData = new Float_t[fNDim];
 123   return;
 124 }
 125
 126 //_________________________________________________________________
 127 void TKDNodeInfo::Bootstrap()
 128 {
 129   if(!fNpar || !fPar) return;
 130
 131   Int_t ndim = Int_t(.5*(TMath::Sqrt(1.+8.*fNpar)-1.))-1;
 132   fNDim = 3*ndim;
 133 }
 134
 135 //_________________________________________________________________
 136 void TKDNodeInfo::SetNode(Int_t ndim, Float_t *data, Float_t *pdf)
 137 {
 138   Build(ndim);
 139   memcpy(fData, data, fNDim*sizeof(Float_t));
 140   fVal[0]=pdf[0]; fVal[1]=pdf[1];
 141 }
 142
 143
 144 //_________________________________________________________________
 145 void TKDNodeInfo::Print(const Option_t *opt) const
 146 {
 147   // Print the content of the node
 148   Int_t dim = Int_t(fNDim/3.);
 149   printf("x[");
 150   for(int idim=0; idim<dim; idim++) printf("%f ", fData?fData[idim]:0.);
 151   printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
 152
 153   if(fData){
 154     Float_t *bounds = &fData[dim];
 155     printf("range[");
 156     for(int idim=0; idim<dim; idim++) printf("(%f %f) ", bounds[2*idim], bounds[2*idim+1]);
 157     printf("]\n");
 158   }
 159   if(strcmp(opt, "a")!=0) return;
 160
 161   if(fNpar){
 162     printf("Fit parameters : \n");
 163     for(int ip=0; ip<fNpar; ip++) printf("p%d[%f] ", ip, fPar[ip]);
 164     printf("\n");
 165   }
 166   if(!fNcov) return;
 167   for(int ip(0), n(0); ip<fNpar; ip++){
 168     for(int jp(ip); jp<fNpar; jp++) printf("c(%d %d)[%f] ", ip, jp, fCov[n++]);
 169     printf("\n");
 170   }
 171 }
 172
 173 //_________________________________________________________________
 174 void TKDNodeInfo::Store(TVectorD const *par, TMatrixD const *cov)
 175 {
 176 // Store the parameters and the covariance in the node
 177
 178   if(!fPar){SetNpar(); fPar = new Double_t[fNpar];}
 179   for(int ip=0; ip<fNpar; ip++) fPar[ip] = (*par)[ip];
 180
 181   if(!cov) return;
 182   if(!fCov){SetNcov(); fCov = new Double_t[fNcov];}
 183   for(int ip(0), np(0); ip<fNpar; ip++)
 184     for(int jp=ip; jp<fNpar; jp++) fCov[np++] = (*cov)(ip,jp);
 185 }
 186
 187 //_________________________________________________________________
 188 Bool_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error) const
 189 {
 190 // Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
 191
 192   Int_t ndim = Int_t(fNDim/3.);
 193   if(ndim>10) return kFALSE; // support only up to 10 dimensions
 194   //printf("ndim[%d] npar[%d] ncov[%d]\n", ndim, fNpar, fNcov);
 195
 196   Double_t fdfdp[66]; memset(fdfdp, 0, ndim*sizeof(Double_t));
 197   Int_t ipar = 0;
 198   fdfdp[ipar++] = 1.;
 199   for(int idim=0; idim<ndim; idim++){
 200     fdfdp[ipar++] = point[idim];
 201     for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
 202   }
 203
 204   // calculate estimation
 205   result =0.; error = 0.;
 206   for(int i=0; i<fNpar; i++) result += fdfdp[i]*fPar[i];
 207   if(!fNcov) return kTRUE;
 208
 209   for(int i(0), n(0); i<fNpar; i++){
 210     error += fdfdp[i]*fdfdp[i]*fCov[n++];
 211     for(int j(i+1); j<fNpar; j++) error += 2.*fdfdp[i]*fdfdp[j]*fCov[n++];
 212   }
 213   error = TMath::Sqrt(error);
 214
 215   //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
 216
 217   return kTRUE;
 218 }
 219
 220
 221
 222 //_________________________________________________________________
 223 TKDNodeInfo::TKDNodeDraw::TKDNodeDraw()
 224   :TBox()
 225   ,fCOG()
 226   ,fNode(NULL)
 227 {
 228   SetFillStyle(3002);
 229   SetFillColor(50+Int_t(gRandom->Uniform()*50.));
 230
 231   fCOG.SetMarkerStyle(3);
 232   fCOG.SetMarkerSize(.7);
 233   fCOG.SetMarkerColor(2);
 234 }
 235
 236
 237 //_________________________________________________________________
 238 void TKDNodeInfo::TKDNodeDraw::Draw(Option_t* option)
 239 {
 240   TBox::Draw(option);
 241   fCOG.Draw("p");
 242 }
 243
 244 //_________________________________________________________________
 245 void TKDNodeInfo::TKDNodeDraw::SetNode(TKDNodeInfo *node, UChar_t size, UChar_t ax1, UChar_t ax2)
 246 {
 247   fNode=node;
 248   const Float_t kBorder = 0.;//1.E-4;
 249   Float_t *bounds = &(node->Data()[size]);
 250   fX1=bounds[2*ax1]+kBorder;
 251   fX2=bounds[2*ax1+1]-kBorder;
 252   fY1=bounds[2*ax2]+kBorder;
 253   fY2=bounds[2*ax2+1]-kBorder;
 254
 255   Float_t x(node->Data()[ax1]), y(node->Data()[ax2]);
 256   fCOG.SetX(x); fCOG.SetY(y);
 257 }
 258
 259
 260 //_________________________________________________________________
 261 void TKDNodeInfo::TKDNodeDraw::Print(const Option_t* option) const
 262 {
 263   if(!fNode) return;
 264   fNode->Print(option);
 265 }