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