PWG2/UNICOR/AliDHN.cxx

   1 // Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2007
   2
   3 //=============================================================================
   4 // multidimensional histogram
   5 // Physically, the data are kept in one single one-dimensional histogram.
   6 // Projecting on 1, 2, and n-1 dimensions is implemented.
   7 // The histogram can be saved on root file as the one-dimensional data
   8 // histogram and the axes, thus eternal forward compatibility is ensured.
   9 //=============================================================================
  10
  11 #include <cmath>
  12 #include <stdlib.h>
  13
  14 #include <TFile.h>
  15 #include <TDirectory.h>
  16 #include <TAxis.h>
  17 #include <TH2.h>
  18 #include "AliDHN.h"
  19
  20 ClassImp(AliDHN)
  21
  22 //=============================================================================
  23 AliDHN::AliDHN(Char_t *nam, Int_t ndim, TAxis **ax)
  24   : TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5),
  25     fNdim(ndim)
  26 {
  27   // Constructor for building from scratch.
  28
  29   // Above, we just have managed to create a single one-dimensional histogram
  30   // with number of bins equal to the product of the numbers of bins in all
  31   // dimensions. For easy inspection the histogram range was set to -0.5,n-0.5.
  32
  33   for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]);
  34
  35   // for speed reasons, number of bins of each axis is stored on fNbins
  36   // and the running product of the latter on fMbins
  37   // so fMbins = {...,fNbins[fNdim-2]*fNbins[fNdim-1],fNbins[fNdim-1],1}
  38
  39   fMbins[fNdim-1] = 1;
  40   for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
  41   for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
  42   printf("%d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX());
  43 }
  44 //=============================================================================
  45 AliDHN::AliDHN(Char_t *filnam, Char_t *nam)
  46   : TH1D(*((TH1D*) TFile::Open(filnam,"read")->GetDirectory(nam)->Get("histo"))),
  47     fNdim(0)
  48 {
  49   // Constructor for reading from file.
  50
  51   TFile *f = TFile::Open(filnam,"read");
  52   f->cd(nam);
  53   TAxis *ax[fMaxNdim];
  54   char axnam[1000];
  55   for (fNdim=0; fNdim<fMaxNdim; fNdim++) {
  56     sprintf(axnam,"axis%d",fNdim);
  57     ax[fNdim] = (TAxis *) gDirectory->Get(axnam);
  58     if (ax[fNdim]) ax[fNdim]->Copy(fAxis[fNdim]);
  59     else break;
  60   }
  61   f->Close();
  62
  63   fMbins[fNdim-1] = 1;
  64   for (int i=0; i<fNdim; i++) fNbins[i] = fAxis[i].GetNbins();
  65   for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
  66
  67   if (GetNbinsX()!=Albins(fNdim,ax)) {
  68     printf("number of bins of histo %d differs from product of nbins of axes %d\n",
  69            GetNbinsX(),Albins(fNdim,ax));
  70     printf("bombing\n");
  71     exit(-1);
  72   }
  73   printf("%d-dimensional histogram %s with %d bins read from file %s\n",
  74          fNdim,nam,GetNbinsX(),filnam);
  75 }
  76 //=============================================================================
  77 Int_t AliDHN::Albins(Int_t n, TAxis **ax)
  78 {
  79   // Calculate product of nbins of ax[0]...ax[n-1]
  80   // (= total number of bins of the multidimensional histogram to be made).
  81
  82   Int_t nbins = 1;
  83   //  while (n--) nbins *= ax[n]->GetNbins();
  84   for (int i=0; i<n; i++) nbins *= ax[i]->GetNbins();
  85   return nbins;
  86 }
  87 //=============================================================================
  88 Int_t AliDHN::MulToOne(Int_t *k) const
  89 {
  90   // Calculate the 1-dim index n from n-dim indices k[fNdim].
  91   // Valid k[i] should be between 0 and fNbins[i]-1.
  92   // Valid result will be between 0 and GetNbinsX()-1.
  93   // Return -1 if under- or overflow in any dimension.
  94
  95   Int_t n = 0;
  96   for (int i=0; i<fNdim; i++) {
  97     if (k[i]<0) return -1;
  98     if (k[i]>=fNbins[i]) return -1;
  99     n += fMbins[i]*k[i];
 100   }
 101   return n;
 102 }
 103 //=============================================================================
 104 Int_t AliDHN::MulToOne(Double_t *x)
 105 {
 106   // Calculate the 1-dim index n from n-dim vector x, representing the
 107   // abscissa of the n-dim histogram. The result will be between 0 and
 108   // GetNbinsX()-1.
 109
 110   Int_t k[fMaxNdim];
 111   for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindBin(x[i])-1;
 112   return MulToOne(k);
 113 }
 114 //=============================================================================
 115 void AliDHN::OneToMul(Int_t n, Int_t *k) const
 116 {
 117   // Calculate the n-dim indices k[fNdim] from 1-dim index n.
 118   // Valid n should be between 0 and GetNbinsX()-1.
 119   // Valid results will be between 0 and fNbins[i]-1.
 120
 121   div_t idi; // integer division structure
 122   for (int i=0; i<fNdim; i++) {
 123     idi  = div(n,fMbins[i]);
 124     k[i] = idi.quot;  // quotient
 125     n    = idi.rem;   // remainder
 126   }
 127 }
 128 //=============================================================================
 129 Int_t AliDHN::Fill(Double_t *xx, Double_t w)
 130 {
 131   // Fill the histogram. The array xx holds the abscissa information, w is the
 132   // weigth. The 1-dim histogram is filled using the standard Fill method
 133   // (direct access to the arrays was tried and was not faster).
 134
 135   int nbin = MulToOne(xx);
 136   if (nbin == -1) return 0;
 137   return TH1D::Fill(nbin+1,w);
 138 }
 139 //=============================================================================
 140 Int_t AliDHN::Fill(Double_t x0, Double_t x1, ...)
 141 {
 142   // Fill the histogram. Arguments are passed as doubles rather than array.
 143   va_list ap;
 144   Double_t xx[fMaxNdim] = {x0, x1};
 145   va_start(ap,x1);
 146   for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t);
 147   Double_t weigth = va_arg(ap,Double_t);
 148   va_end(ap);
 149   return Fill(xx,weigth);
 150 }
 151 //=============================================================================
 152 Int_t AliDHN::Write() const
 153 {
 154   // Save the 1-dim histo and the axes in a subdirectory on file. This might
 155   // not be the most elegant way but it is very simple and backward and forward
 156   // compatible.
 157
 158   Int_t nbytes = 0;
 159   TH1D histo(*this);
 160   histo.SetName("histo"); // hadd bug fix
 161   TDirectory *dest = gDirectory->mkdir(GetName());
 162   if (dest) {
 163     dest->cd();
 164     nbytes += histo.Write("histo");
 165     char axnam[1000];
 166     for (int i=0; i<fNdim; i++) {
 167       sprintf(axnam,"axis%d",i);
 168       nbytes += fAxis[i].Write(axnam);
 169     }
 170     printf("   %s stored in %s\n",GetName(),dest->GetPath());
 171     dest->cd("..");
 172   }
 173   return nbytes;
 174 }
 175 //=============================================================================
 176 AliDHN *AliDHN::ProjectAlong(char *nam, Int_t dim, Int_t first, Int_t last)
 177 {
 178   // Reduce dimension dim by summing up its bins between first and last.
 179   // Use root convention: bin=1 is the first bin, bin=nbins is the last.
 180   // Return the resulting fNdim-1 dimensional histogram.
 181
 182   if (dim<0 || dim>fNdim-1) return 0;
 183   if (first<0) first = 1;
 184   if (last<0) last = fNbins[dim];
 185
 186   // create new (reduced) histogram
 187
 188   TAxis *ax[fMaxNdim-1];
 189   int n=0;
 190   for (int i=0; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i);
 191   char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim);
 192   char *eame = Form("%s_error",name);
 193   AliDHN *his = new AliDHN(name,fNdim-1,ax); // result histogram
 194   AliDHN *eis = new AliDHN(eame,fNdim-1,ax); // temporary storage for errors
 195
 196   // sum up the content and errors squared
 197
 198   Int_t *k = new Int_t[fNdim];   // old hist multiindex
 199   Int_t *m = new Int_t[fNdim-1]; // new hist multiindex
 200   for (int i=0; i<GetNbinsX(); i++) {
 201     OneToMul(i,k);
 202     if (k[dim]+1<first) continue;
 203     if (k[dim]+1>last) continue;
 204     n = 0;
 205     for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j];
 206     n = his->MulToOne(m);
 207     his->AddBinContent(n+1,GetBinContent(i+1));
 208     eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1));
 209   }
 210
 211   // combine content and errors in one histogram
 212
 213   for (int i=0; i<his->GetNbinsX(); i++) {
 214     his->SetBinError(i+1,sqrt(eis->GetBinContent(i+1)));
 215   }
 216
 217   his->SetLineColor(this->GetLineColor());
 218   his->SetFillColor(this->GetFillColor());
 219   his->SetMarkerColor(this->GetMarkerColor());
 220   his->SetMarkerStyle(this->GetMarkerStyle());
 221
 222   // some cleanup
 223
 224   delete eis;
 225   return his;
 226 }
 227 //=============================================================================
 228 TH1D *AliDHN::ProjectOn(char *nam, Int_t dim, Int_t *first, Int_t *last) const
 229 {
 230   // Project on dimension dim. Use only bins between first[i] and last[i].
 231   // Use root convention: bin=1 is the first bin, bin=nbins is the last.
 232   // first[i]=0 means from the first bin
 233   // last[i]=0 means till the last bin
 234
 235   if (dim<0 || dim>fNdim-1) return 0;
 236
 237   // calculate the projection; lowest index 0
 238
 239   double *yy = new double[fNbins[dim]]; // value
 240   double *ey = new double[fNbins[dim]]; // error
 241   for (int i=0; i<fNbins[dim]; i++) yy[i]=0;
 242   for (int i=0; i<fNbins[dim]; i++) ey[i]=0;
 243   Int_t *k = new Int_t[fNdim];
 244   for (int i=0; i<GetNbinsX(); i++) {
 245     OneToMul(i,k);
 246     int isgood = 1; // bin within the range?
 247     for (int j=0; j<fNdim; j++) {
 248       if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
 249       if (last)  if (last[j])  if (k[j]+1>last[j])  {isgood = 0; break;}
 250     }
 251     if (isgood) {
 252       yy[k[dim]]+=GetBinContent(i+1);
 253       ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1);
 254     }
 255   }
 256
 257   // make the projection histogram
 258   // use name nam if specified; otherwise generate one
 259
 260   TH1D *his;
 261   char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
 262   if (fAxis[dim].IsVariableBinSize())
 263     his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray());
 264   else
 265     his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax());
 266   his->SetXTitle(fAxis[dim].GetTitle());
 267   // or his->GetXaxis()->ImportAttributes(ax);
 268   his->Sumw2();
 269   his->SetLineColor(this->GetLineColor());
 270   his->SetFillColor(this->GetFillColor());
 271   his->SetMarkerColor(this->GetMarkerColor());
 272   his->SetMarkerStyle(this->GetMarkerStyle());
 273   for (int i=0; i<his->GetNbinsX(); i++) {
 274     his->SetBinContent(i+1,yy[i]);
 275     his->SetBinError(i+1,sqrt(ey[i]));
 276   }
 277
 278   // some cleanup
 279
 280   delete [] yy;
 281   delete [] ey;
 282   delete [] k;
 283   //  if (name!=nam) delete [] name;
 284
 285   return his;
 286 }
 287 //=============================================================================
 288 TH1D *AliDHN::ProjectOn(char *nam, Int_t dim, Double_t *first, Double_t *last)
 289 {
 290   // Project on dimension dim. Use only bins between first[i] and last[i].
 291
 292   if (dim<0 || dim>fNdim-1) return 0;
 293   Int_t kfirst[fMaxNdim];
 294   Int_t klast[fMaxNdim];
 295   for (int i=0; i<fNdim; i++) {
 296     kfirst[i] = fAxis[i].FindBin(first[i]);
 297     klast[i] = fAxis[i].FindBin(last[i]);
 298   }
 299   return ProjectOn(nam,dim,kfirst,klast);
 300 }
 301 //=============================================================================
 302 TH2D *AliDHN::ProjectOn(char *nam, Int_t dim0, Int_t dim1, Int_t *first, Int_t *last) const
 303 {
 304   // Project on dim1 vs dim0. Use only bins between first[i] and last[i].
 305   // Use root convention: bin=1 is the first bin, bin=nbins is the last.
 306   // first[i]=0 means from the first bin
 307   // last[i]=0 means till the last bin
 308
 309   if (dim0<0 || dim0>fNdim-1) return 0;
 310   if (dim1<0 || dim1>fNdim-1) return 0;
 311
 312   // calculate the projection
 313
 314   double **yy = new double*[fNbins[dim0]]; // value
 315   double **ey = new double*[fNbins[dim0]]; // error
 316   for (int i=0; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]];
 317   for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]];
 318   for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0;
 319   for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0;
 320   Int_t *k = new Int_t[fNdim];
 321   for (int i=0; i<GetNbinsX(); i++) {
 322     OneToMul(i,k);
 323     int isgood = 1; // bin within the range?
 324     for (int j=0; j<fNdim; j++) {
 325       if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
 326       if (last)  if (last[j])  if (k[j]+1>last[j])  {isgood = 0; break;}
 327     }
 328     if (isgood) {
 329       yy[k[dim0]][k[dim1]]+=GetBinContent(i+1);
 330       ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1);
 331     }
 332   }
 333
 334   // make the projection histogram
 335   // use name nam if specified; otherwise generate one
 336
 337   TH2D *his=0;
 338   char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0);
 339   if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
 340     his = new TH2D(name,name,
 341                    fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
 342                    fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
 343   else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
 344     his = new TH2D(name,name,
 345                    fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
 346                    fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
 347   else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
 348     his = new TH2D(name,name,
 349                    fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
 350                    fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
 351   else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
 352     his = new TH2D(name,name,
 353                    fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
 354                    fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
 355   his->SetXTitle(fAxis[dim0].GetTitle());
 356   his->SetYTitle(fAxis[dim1].GetTitle());
 357   // or his->GetXaxis()->ImportAttributes(ax0);
 358   // or his->GetYaxis()->ImportAttributes(ax1);
 359   his->Sumw2();
 360   his->SetLineColor(this->GetLineColor());
 361   his->SetFillColor(this->GetFillColor());
 362   his->SetMarkerColor(this->GetMarkerColor());
 363   his->SetMarkerStyle(this->GetMarkerStyle());
 364   for (int i=0; i<his->GetNbinsX(); i++)
 365   for (int j=0; j<his->GetNbinsY(); j++) {
 366     his->SetBinContent(i+1,j+1,yy[i][j]);
 367     his->SetBinError(i+1,j+1,sqrt(ey[i][j]));
 368   }
 369
 370   // some cleanup
 371
 372   for (int i=0; i<fNbins[dim0]; i++) delete [] yy[i];
 373   for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i];
 374   delete [] yy;
 375   delete [] ey;
 376   delete [] k;
 377   //  if (name!=nam) delete [] name;
 378
 379   return his;
 380 }
 381 //=============================================================================