#include "AliTHn.h"
#include "TList.h"
#include "TCollection.h"
#include "AliLog.h"
#include "TArrayF.h"
#include "THnSparse.h"
#include "TMath.h"
ClassImp(AliTHn)
AliTHn::AliTHn() :
AliCFContainer(),
fNBins(0),
fNVars(0),
fNSteps(0),
fValues(0),
fSumw2(0),
axisCache(0),
fNbinsCache(0),
fLastVars(0),
fLastBins(0)
{
}
AliTHn::AliTHn(const Char_t* name, const Char_t* title,const Int_t nSelStep, const Int_t nVarIn, const Int_t* nBinIn) :
AliCFContainer(name, title, nSelStep, nVarIn, nBinIn),
fNBins(0),
fNVars(nVarIn),
fNSteps(nSelStep),
fValues(0),
fSumw2(0),
axisCache(0),
fNbinsCache(0),
fLastVars(0),
fLastBins(0)
{
fNBins = 1;
for (Int_t i=0; i<fNVars; i++)
fNBins *= nBinIn[i];
Init();
}
void AliTHn::Init()
{
fValues = new TArrayF*[fNSteps];
fSumw2 = new TArrayF*[fNSteps];
for (Int_t i=0; i<fNSteps; i++)
{
fValues[i] = 0;
fSumw2[i] = 0;
}
}
AliTHn::AliTHn(const AliTHn &c) :
AliCFContainer(c),
fNBins(c.fNBins),
fNVars(c.fNVars),
fNSteps(c.fNSteps),
fValues(new TArrayF*[c.fNSteps]),
fSumw2(new TArrayF*[c.fNSteps]),
axisCache(0),
fNbinsCache(0),
fLastVars(0),
fLastBins(0)
{
memset(fValues,0,fNSteps*sizeof(TArrayF*));
memset(fSumw2,0,fNSteps*sizeof(TArrayF*));
for (Int_t i=0; i<fNSteps; i++) {
if (c.fValues[i]) fValues[i] = new TArrayF(*(c.fValues[i]));
if (c.fSumw2[i]) fSumw2[i] = new TArrayF(*(c.fSumw2[i]));
}
}
AliTHn::~AliTHn()
{
DeleteContainers();
delete[] fValues;
delete[] fSumw2;
delete[] axisCache;
delete[] fNbinsCache;
delete[] fLastVars;
delete[] fLastBins;
}
void AliTHn::DeleteContainers()
{
for (Int_t i=0; i<fNSteps; i++)
{
if (fValues && fValues[i])
{
delete fValues[i];
fValues[i] = 0;
}
if (fSumw2 && fSumw2[i])
{
delete fSumw2[i];
fSumw2[i] = 0;
}
}
}
AliTHn &AliTHn::operator=(const AliTHn &c)
{
if (this != &c) {
AliCFContainer::operator=(c);
fNBins=c.fNBins;
fNVars=c.fNVars;
if(fNSteps) {
for(Int_t i=0; i< fNSteps; ++i) {
delete fValues[i];
delete fSumw2[i];
}
delete [] fValues;
delete [] fSumw2;
}
fNSteps=c.fNSteps;
if(fNSteps) {
fValues=new TArrayF*[fNSteps];
fSumw2=new TArrayF*[fNSteps];
memset(fValues,0,fNSteps*sizeof(TArrayF*));
memset(fSumw2,0,fNSteps*sizeof(TArrayF*));
for (Int_t i=0; i<fNSteps; i++) {
if (c.fValues[i]) fValues[i] = new TArrayF(*(c.fValues[i]));
if (c.fSumw2[i]) fSumw2[i] = new TArrayF(*(c.fSumw2[i]));
}
} else {
fValues = 0;
fSumw2 = 0;
}
delete [] axisCache;
axisCache = new TAxis*[fNVars];
memcpy(axisCache, c.axisCache, fNVars*sizeof(TAxis*));
}
return *this;
}
void AliTHn::Copy(TObject& c) const
{
AliTHn& target = (AliTHn &) c;
AliCFContainer::Copy(target);
target.fNSteps = fNSteps;
target.fNBins = fNBins;
target.fNVars = fNVars;
target.Init();
for (Int_t i=0; i<fNSteps; i++)
{
if (fValues[i])
target.fValues[i] = new TArrayF(*(fValues[i]));
else
target.fValues[i] = 0;
if (fSumw2[i])
target.fSumw2[i] = new TArrayF(*(fSumw2[i]));
else
target.fSumw2[i] = 0;
}
}
Long64_t AliTHn::Merge(TCollection* list)
{
if (!list)
return 0;
if (list->IsEmpty())
return 1;
AliCFContainer::Merge(list);
TIterator* iter = list->MakeIterator();
TObject* obj;
Int_t count = 0;
while ((obj = iter->Next())) {
AliTHn* entry = dynamic_cast<AliTHn*> (obj);
if (entry == 0)
continue;
for (Int_t i=0; i<fNSteps; i++)
{
if (entry->fValues[i])
{
if (!fValues[i])
fValues[i] = new TArrayF(fNBins);
for (Long64_t l = 0; l<fNBins; l++)
fValues[i]->GetArray()[l] += entry->fValues[i]->GetArray()[l];
}
if (entry->fSumw2[i])
{
if (!fSumw2[i])
fSumw2[i] = new TArrayF(fNBins);
for (Long64_t l = 0; l<fNBins; l++)
fSumw2[i]->GetArray()[l] += entry->fSumw2[i]->GetArray()[l];
}
}
count++;
}
return count+1;
}
void AliTHn::Fill(const Double_t *var, Int_t istep, Double_t weight)
{
if (!axisCache)
{
axisCache = new TAxis*[fNVars];
fNbinsCache = new Int_t[fNVars];
for (Int_t i=0; i<fNVars; i++)
{
axisCache[i] = GetAxis(i, 0);
fNbinsCache[i] = axisCache[i]->GetNbins();
}
fLastVars = new Double_t[fNVars];
fLastBins = new Int_t[fNVars];
for (Int_t i=0; i<fNVars; i++)
{
fLastBins[i] = axisCache[i]->FindBin(var[i]);
fLastVars[i] = var[i];
}
}
Long64_t bin = 0;
for (Int_t i=0; i<fNVars; i++)
{
bin *= fNbinsCache[i];
Int_t tmpBin = 0;
if (fLastVars[i] == var[i])
tmpBin = fLastBins[i];
else
{
tmpBin = axisCache[i]->FindBin(var[i]);
fLastBins[i] = tmpBin;
fLastVars[i] = var[i];
}
if (tmpBin < 1 || tmpBin > fNbinsCache[i])
return;
bin += tmpBin - 1;
}
if (!fValues[istep])
{
fValues[istep] = new TArrayF(fNBins);
AliInfo(Form("Created values container for step %d", istep));
}
if (weight != 1)
{
if (!fSumw2[istep])
{
fSumw2[istep] = new TArrayF(*fValues[istep]);
AliInfo(Form("Created sumw2 container for step %d", istep));
}
}
fValues[istep]->GetArray()[bin] += weight;
if (fSumw2[istep])
fSumw2[istep]->GetArray()[bin] += weight * weight;
}
Long64_t AliTHn::GetGlobalBinIndex(const Int_t* binIdx)
{
Long64_t bin = 0;
for (Int_t i=0; i<fNVars; i++)
{
bin *= GetAxis(i, 0)->GetNbins();
bin += binIdx[i] - 1;
}
return bin;
}
void AliTHn::FillContainer(AliCFContainer* cont)
{
for (Int_t i=0; i<fNSteps; i++)
{
if (!fValues[i])
continue;
Float_t* source = fValues[i]->GetArray();
Float_t* sourceSumw2 = source;
if (fSumw2[i])
sourceSumw2 = fSumw2[i]->GetArray();
THnSparse* target = cont->GetGrid(i)->GetGrid();
Int_t* binIdx = new Int_t[fNVars];
Int_t* nBins = new Int_t[fNVars];
for (Int_t j=0; j<fNVars; j++)
{
binIdx[j] = 1;
nBins[j] = target->GetAxis(j)->GetNbins();
}
Long64_t count = 0;
while (1)
{
Long64_t globalBin = GetGlobalBinIndex(binIdx);
if (source[globalBin] != 0)
{
target->SetBinContent(binIdx, source[globalBin]);
target->SetBinError(binIdx, TMath::Sqrt(sourceSumw2[globalBin]));
count++;
}
binIdx[fNVars-1]++;
for (Int_t j=fNVars-1; j>0; j--)
{
if (binIdx[j] > nBins[j])
{
binIdx[j] = 1;
binIdx[j-1]++;
}
}
if (binIdx[0] > nBins[0])
break;
}
AliInfo(Form("Step %d: copied %lld entries out of %lld bins", i, count, GetGlobalBinIndex(binIdx)));
delete[] binIdx;
delete[] nBins;
}
}
void AliTHn::FillParent()
{
FillContainer(this);
}
void AliTHn::ReduceAxis()
{
Int_t axis = fNVars-1;
for (Int_t i=0; i<fNSteps; i++)
{
if (!fValues[i])
continue;
Float_t* source = fValues[i]->GetArray();
Float_t* sourceSumw2 = 0;
if (fSumw2[i])
sourceSumw2 = fSumw2[i]->GetArray();
THnSparse* target = GetGrid(i)->GetGrid();
Int_t* binIdx = new Int_t[fNVars];
Int_t* nBins = new Int_t[fNVars];
for (Int_t j=0; j<fNVars; j++)
{
binIdx[j] = 1;
nBins[j] = target->GetAxis(j)->GetNbins();
}
Long64_t count = 0;
while (1)
{
Float_t sumValues = 0;
Float_t sumSumw2 = 0;
for (Int_t j=1; j<=nBins[axis]; j++)
{
binIdx[axis] = j;
Long64_t globalBin = GetGlobalBinIndex(binIdx);
sumValues += source[globalBin];
source[globalBin] = 0;
if (sourceSumw2)
{
sumSumw2 += sourceSumw2[globalBin];
sourceSumw2[globalBin] = 0;
}
}
binIdx[axis] = 1;
Long64_t globalBin = GetGlobalBinIndex(binIdx);
source[globalBin] = sumValues;
if (sourceSumw2)
sourceSumw2[globalBin] = sumSumw2;
count++;
binIdx[fNVars-2]++;
for (Int_t j=fNVars-2; j>0; j--)
{
if (binIdx[j] > nBins[j])
{
binIdx[j] = 1;
binIdx[j-1]++;
}
}
if (binIdx[0] > nBins[0])
break;
}
AliInfo(Form("Step %d: reduced %lld bins to %lld entries", i, GetGlobalBinIndex(binIdx), count));
delete[] binIdx;
delete[] nBins;
}
}