#include <cmath>
#include <stdlib.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TROOT.h>
#include <TAxis.h>
#include <TH2.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TString.h>
#include "AliUnicorHN.h"
ClassImp(AliUnicorHN)
AliUnicorHN::AliUnicorHN(const char *nam, Int_t ndim, TAxis **ax)
: TH1D(nam, nam, Albins(ndim,ax), 0.5, Albins(ndim,ax)+0.5),
fNdim(ndim) {
for (int i=0; i<fNdim; i++) ax[i]->Copy(fAxis[i]);
for (int i=0; i<fNdim; i++) fAxis[i].SetName(Form("axis%d",i));
for (int i=0; i<fgkMaxNdim; i++) fNbins[i] = fAxis[i].GetNbins();
for (int i=0; i<fgkMaxNdim; i++) fMbins[i] = 1;
for (int i=fNdim-1; i>0; i--) fMbins[i-1] = fMbins[i]*fNbins[i];
printf(" %d-dimensional histogram %s with %d bins created\n",fNdim,nam,GetNbinsX());
}
AliUnicorHN* AliUnicorHN::Retrieve(const char *filnam, const char *nam) {
TFile *f = TFile::Open(filnam,"read");
if (!f) return 0;
if (!f->cd(nam)) {f->Close(); return 0;}
TH1D *hist = (TH1D*) gDirectory->Get("histo");
if (!hist) {printf("No histogram histo on file %s in directory %s\n",filnam,nam); return 0;}
hist->SetDirectory(gROOT);
TAxis *ax[fgkMaxNdim]={0};
int n=0;
while ((ax[n] = (TAxis *) gDirectory->Get(Form("axis%d",n)))) n++;
f->Close();
if (hist->GetNbinsX()!=Albins(n,ax)) {
printf("number of bins of histo %d differs from product of nbins of axes %d\n",
hist->GetNbinsX(),Albins(n,ax));
return 0;
}
TString path=nam;
TObjArray *ar = path.Tokenize("/");
TObjString *os = (TObjString*) ar->Last();
const char *lastnam = os->GetString().Data();
AliUnicorHN *hi = new AliUnicorHN(lastnam,n,ax);
hist->Copy(*((TH1D*)hi));
hi->SetName(lastnam);
delete hist;
delete ar;
return hi;
}
Int_t AliUnicorHN::Albins(Int_t n, TAxis **ax) {
Int_t nbins = 1;
for (int i=0; i<n; i++) nbins *= ax[i]->GetNbins();
return nbins;
}
Int_t AliUnicorHN::MulToOne(const Int_t * const k) const {
Int_t n = 0;
for (int i=0; i<fNdim; i++) {
if (k[i]<0) return -1;
if (k[i]>=fNbins[i]) return -1;
n += fMbins[i]*k[i];
}
return n;
}
Int_t AliUnicorHN::MulToOne(Double_t *x) {
Int_t k[fgkMaxNdim] = {0};
for (int i=0; i<fNdim; i++) k[i] = fAxis[i].FindFixBin(x[i])-1;
return MulToOne(k);
}
void AliUnicorHN::OneToMul(Int_t n, Int_t *k) const {
div_t idi;
for (int i=0; i<fNdim; i++) {
idi = div(n,fMbins[i]);
k[i] = idi.quot;
n = idi.rem;
}
}
Int_t AliUnicorHN::Fill(Double_t *xx, Double_t w) {
int nbin = MulToOne(xx);
if (nbin == -1) return 0;
return TH1D::Fill(nbin+1,w);
}
Int_t AliUnicorHN::Fill(Double_t x0, Double_t x1, ...) {
va_list ap;
Double_t xx[fgkMaxNdim] = {x0, x1};
va_start(ap,x1);
for (int i=2; i<fNdim; i++) xx[i] = va_arg(ap,Double_t);
Double_t weigth = va_arg(ap,Double_t);
va_end(ap);
return Fill(xx,weigth);
}
Int_t AliUnicorHN::Save() const {
Int_t nbytes = 0;
TH1D histo(*this);
histo.SetName("histo");
TDirectory *dest = gDirectory->mkdir(GetName());
if (dest) {
dest->cd();
nbytes += histo.Write("histo");
for (int i=0; i<fNdim; i++) nbytes += fAxis[i].Write();
printf(" %s stored in %s\n",GetName(),dest->GetPath());
dest->cd("..");
}
return nbytes;
}
AliUnicorHN *AliUnicorHN::ProjectAlong(const char *nam, Int_t dim, Int_t first, Int_t last) {
if (dim<0 || dim>fNdim-1) return 0;
if (last<=0) last = fNbins[dim];
TAxis *ax[fgkMaxNdim-1] = {0};
int n=0;
for (int i=0; i<fNdim; i++) if (i!=dim) ax[n++] = GetAxis(i);
const char *name = strlen(nam)? nam : Form("%s_wo%d",GetName(),dim);
const char *eame = Form("%s_error",name);
AliUnicorHN *his = new AliUnicorHN(name,fNdim-1,ax);
AliUnicorHN *eis = new AliUnicorHN(eame,fNdim-1,ax);
int k[fgkMaxNdim] = {0};
int m[fgkMaxNdim] = {0};
for (int i=0; i<GetNbinsX(); i++) {
OneToMul(i,k);
if (k[dim]+1<first) continue;
if (k[dim]+1>last) continue;
n = 0;
for (int j=0; j<fNdim; j++) if (j!=dim) m[n++] = k[j];
n = his->MulToOne(m);
his->AddBinContent(n+1,GetBinContent(i+1));
eis->AddBinContent(n+1,GetBinError(i+1)*GetBinError(i+1));
}
for (int i=0; i<his->GetNbinsX(); i++) {
his->SetBinError(i+1,sqrt(eis->GetBinContent(i+1)));
}
his->SetLineColor(this->GetLineColor());
his->SetFillColor(this->GetFillColor());
his->SetMarkerColor(this->GetMarkerColor());
his->SetMarkerStyle(this->GetMarkerStyle());
delete eis;
return his;
}
TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Int_t * const first,
const Int_t * const last) const {
if (dim<0 || dim>fNdim-1) return 0;
double *yy = new double[fNbins[dim]];
double *ey = new double[fNbins[dim]];
for (int i=0; i<fNbins[dim]; i++) yy[i]=0;
for (int i=0; i<fNbins[dim]; i++) ey[i]=0;
Int_t *k = new Int_t[fNdim];
for (int i=0; i<GetNbinsX(); i++) {
OneToMul(i,k);
int isgood = 1;
for (int j=0; j<fNdim; j++) {
if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;}
}
if (isgood) {
yy[k[dim]]+=GetBinContent(i+1);
ey[k[dim]]+=GetBinError(i+1)*GetBinError(i+1);
}
}
TH1D *his;
const char *name = strlen(nam)? nam : Form("%s_proj%d",GetName(),dim);
if (fAxis[dim].IsVariableBinSize())
his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXbins()->GetArray());
else
his = new TH1D(name,name,fNbins[dim],fAxis[dim].GetXmin(),fAxis[dim].GetXmax());
his->SetXTitle(fAxis[dim].GetTitle());
his->Sumw2();
his->SetLineColor(this->GetLineColor());
his->SetFillColor(this->GetFillColor());
his->SetMarkerColor(this->GetMarkerColor());
his->SetMarkerStyle(this->GetMarkerStyle());
for (int i=0; i<his->GetNbinsX(); i++) {
his->SetBinContent(i+1,yy[i]);
his->SetBinError(i+1,sqrt(ey[i]));
}
delete [] yy;
delete [] ey;
delete [] k;
return his;
}
TH1D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim, const Double_t * const first,
const Double_t * const last) const {
if (dim<0 || dim>fNdim-1) return 0;
Int_t kfirst[fgkMaxNdim] = {0};
Int_t klast[fgkMaxNdim] = {0};
for (int i=0; i<fNdim; i++) {
kfirst[i] = fAxis[i].FindFixBin(first[i]);
klast[i] = fAxis[i].FindFixBin(last[i]);
}
return ProjectOn(nam,dim,kfirst,klast);
}
TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Int_t *
const first, const Int_t * const last) const {
if (dim0<0 || dim0>fNdim-1) return 0;
if (dim1<0 || dim1>fNdim-1) return 0;
double **yy = new double*[fNbins[dim0]];
double **ey = new double*[fNbins[dim0]];
for (int i=0; i<fNbins[dim0]; i++) yy[i] = new double[fNbins[dim1]];
for (int i=0; i<fNbins[dim0]; i++) ey[i] = new double[fNbins[dim1]];
for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) yy[i][j]=0;
for (int i=0; i<fNbins[dim0]; i++) for (int j=0; j<fNbins[dim1]; j++) ey[i][j]=0;
Int_t *k = new Int_t[fNdim];
for (int i=0; i<GetNbinsX(); i++) {
OneToMul(i,k);
int isgood = 1;
for (int j=0; j<fNdim; j++) {
if (first) if (first[j]) if (k[j]+1<first[j]) {isgood = 0; break;}
if (last) if (last[j]) if (k[j]+1>last[j]) {isgood = 0; break;}
}
if (isgood) {
yy[k[dim0]][k[dim1]]+=GetBinContent(i+1);
ey[k[dim0]][k[dim1]]+=GetBinError(i+1)*GetBinError(i+1);
}
}
TH2D *his=0;
const char *name = strlen(nam)? nam : Form("%s_proj%dvs%d",GetName(),dim1,dim0);
if (fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
his = new TH2D(name,name,
fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
else if (!fAxis[dim0].IsVariableBinSize() && fAxis[dim1].IsVariableBinSize())
his = new TH2D(name,name,
fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
fNbins[dim1],fAxis[dim1].GetXbins()->GetArray());
else if (fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
his = new TH2D(name,name,
fNbins[dim0],fAxis[dim0].GetXbins()->GetArray(),
fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
else if (!fAxis[dim0].IsVariableBinSize() && !fAxis[dim1].IsVariableBinSize())
his = new TH2D(name,name,
fNbins[dim0],fAxis[dim0].GetXmin(),fAxis[dim0].GetXmax(),
fNbins[dim1],fAxis[dim1].GetXmin(),fAxis[dim1].GetXmax());
his->SetXTitle(fAxis[dim0].GetTitle());
his->SetYTitle(fAxis[dim1].GetTitle());
his->Sumw2();
his->SetLineColor(this->GetLineColor());
his->SetFillColor(this->GetFillColor());
his->SetMarkerColor(this->GetMarkerColor());
his->SetMarkerStyle(this->GetMarkerStyle());
for (int i=0; i<his->GetNbinsX(); i++)
for (int j=0; j<his->GetNbinsY(); j++) {
his->SetBinContent(i+1,j+1,yy[i][j]);
his->SetBinError(i+1,j+1,sqrt(ey[i][j]));
}
for (int i=0; i<fNbins[dim0]; i++) delete [] yy[i];
for (int i=0; i<fNbins[dim0]; i++) delete [] ey[i];
delete [] yy;
delete [] ey;
delete [] k;
return his;
}
TH2D *AliUnicorHN::ProjectOn(const char *nam, Int_t dim0, Int_t dim1, const Double_t *
const first, const Double_t * const last) const {
if (dim0<0 || dim0>fNdim-1) return 0;
if (dim1<0 || dim1>fNdim-1) return 0;
Int_t kfirst[fgkMaxNdim] = {0};
Int_t klast[fgkMaxNdim] = {0};
for (int i=0; i<fNdim; i++) {
kfirst[i] = fAxis[i].FindFixBin(first[i]);
klast[i] = fAxis[i].FindFixBin(last[i]);
}
return ProjectOn(nam,dim0,dim1,kfirst,klast);
}