PWG2/UNICOR/AliUnicorHN.cxx

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