#include "AliCFGridSparse.h"
#include "THnSparse.h"
#include "AliLog.h"
#include "TMath.h"
#include "TROOT.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TAxis.h"
#include "AliCFUnfolding.h"
ClassImp(AliCFGridSparse)
AliCFGridSparse::AliCFGridSparse() :
AliCFFrame(),
fSumW2(kFALSE),
fData(0x0)
{
}
AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title) :
AliCFFrame(name,title),
fSumW2(kFALSE),
fData(0x0)
{
}
AliCFGridSparse::AliCFGridSparse(const Char_t* name, const Char_t* title, Int_t nVarIn, const Int_t * nBinIn) :
AliCFFrame(name,title),
fSumW2(kFALSE),
fData(0x0)
{
fData = new THnSparseF(name,title,nVarIn,nBinIn);
}
AliCFGridSparse::~AliCFGridSparse()
{
if (fData) delete fData;
}
AliCFGridSparse::AliCFGridSparse(const AliCFGridSparse& c) :
AliCFFrame(c),
fSumW2(kFALSE),
fData(0x0)
{
((AliCFGridSparse &)c).Copy(*this);
}
AliCFGridSparse& AliCFGridSparse::operator=(const AliCFGridSparse &c)
{
if (this != &c) c.Copy(*this);
return *this;
}
void AliCFGridSparse::SetBinLimits(Int_t ivar, Double_t min, Double_t max)
{
Int_t nBins = GetNBins(ivar);
Double_t * array = new Double_t[nBins+1];
for (Int_t iEdge=0; iEdge<=nBins; iEdge++) array[iEdge] = min + iEdge * (max-min)/nBins ;
fData->SetBinEdges(ivar, array);
delete [] array ;
}
void AliCFGridSparse::SetBinLimits(Int_t ivar, const Double_t *array)
{
fData->SetBinEdges(ivar, array);
}
void AliCFGridSparse::Fill(const Double_t *var, Double_t weight)
{
fData->Fill(var,weight);
}
AliCFGridSparse* AliCFGridSparse::MakeSlice(Int_t nVars, const Int_t* vars, const Double_t* varMin, const Double_t* varMax, Bool_t useBins) const
{
Int_t* bins = new Int_t[nVars];
for (Int_t iVar=0; iVar<nVars; iVar++) {
bins[iVar] = GetNBins(vars[iVar]);
}
AliCFGridSparse* out = new AliCFGridSparse(fName,fTitle,nVars,bins);
THnSparse* clone = ((THnSparse*)fData->Clone());
if (varMin && varMax) {
for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) {
SetAxisRange(clone->GetAxis(iAxis),varMin[iAxis],varMax[iAxis],useBins);
}
}
else AliInfo("Keeping same axis ranges");
out->SetGrid(clone->Projection(nVars,vars));
delete [] bins;
delete clone;
return out;
}
Float_t AliCFGridSparse::GetBinCenter(Int_t ivar, Int_t ibin) const
{
return (Float_t) fData->GetAxis(ivar)->GetBinCenter(ibin);
}
Float_t AliCFGridSparse::GetBinSize(Int_t ivar, Int_t ibin) const
{
return (Float_t) fData->GetAxis(ivar)->GetBinUpEdge(ibin) - fData->GetAxis(ivar)->GetBinLowEdge(ibin);
}
Float_t AliCFGridSparse::GetEntries() const
{
return fData->GetEntries();
}
Float_t AliCFGridSparse::GetElement(Long_t index) const
{
return fData->GetBinContent(index);
}
Float_t AliCFGridSparse::GetElement(const Int_t *bin) const
{
return fData->GetBinContent(bin);
}
Float_t AliCFGridSparse::GetElement(const Double_t *var) const
{
Long_t index = fData->GetBin(var,kFALSE);
if (index<0) return 0.;
return fData->GetBinContent(index);
}
Float_t AliCFGridSparse::GetElementError(Long_t index) const
{
return fData->GetBinError(index);
}
Float_t AliCFGridSparse::GetElementError(const Int_t *bin) const
{
return fData->GetBinError(bin);
}
Float_t AliCFGridSparse::GetElementError(const Double_t *var) const
{
Long_t index=fData->GetBin(var,kFALSE);
if (index<0) return 0.;
return fData->GetBinError(index);
}
void AliCFGridSparse::SetElement(Long_t index, Float_t val)
{
Int_t* bin = new Int_t[GetNVar()];
fData->GetBinContent(index,bin);
SetElement(bin,val);
delete [] bin ;
}
void AliCFGridSparse::SetElement(const Int_t *bin, Float_t val)
{
fData->SetBinContent(bin,val);
}
void AliCFGridSparse::SetElement(const Double_t *var, Float_t val)
{
Long_t index=fData->GetBin(var,kTRUE);
Int_t *bin = new Int_t[GetNVar()];
fData->GetBinContent(index,bin);
SetElement(bin,val);
delete [] bin;
}
void AliCFGridSparse::SetElementError(Long_t index, Float_t val)
{
Int_t *bin = new Int_t[GetNVar()];
fData->GetBinContent(index,bin);
SetElementError(bin,val);
delete [] bin;
}
void AliCFGridSparse::SetElementError(const Int_t *bin, Float_t val)
{
fData->SetBinError(bin,val);
}
void AliCFGridSparse::SetElementError(const Double_t *var, Float_t val)
{
Long_t index=fData->GetBin(var);
Int_t *bin = new Int_t[GetNVar()];
fData->GetBinContent(index,bin);
SetElementError(bin,val);
delete [] bin;
}
void AliCFGridSparse::SumW2()
{
if(!fSumW2){
fData->CalculateErrors(kTRUE);
}
fSumW2=kTRUE;
}
void AliCFGridSparse::Add(const AliCFGridSparse* aGrid, Double_t c)
{
if (aGrid->GetNVar() != GetNVar()){
AliError("Different number of variables, cannot add the grids");
return;
}
if (!fSumW2 && aGrid->GetSumW2()) SumW2();
fData->Add(aGrid->GetGrid(),c);
}
void AliCFGridSparse::Add(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
{
if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
AliInfo("Different number of variables, cannot add the grids");
return;
}
if (!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
fData->Reset();
fData->Add(aGrid1->GetGrid(),c1);
fData->Add(aGrid2->GetGrid(),c2);
}
void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid, Double_t c)
{
if (aGrid->GetNVar() != GetNVar()) {
AliError("Different number of variables, cannot multiply the grids");
return;
}
if(!fSumW2 && aGrid->GetSumW2()) SumW2();
THnSparse *h = aGrid->GetGrid();
fData->Multiply(h);
fData->Scale(c);
}
void AliCFGridSparse::Multiply(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2)
{
if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
AliError("Different number of variables, cannot multiply the grids");
return;
}
if(!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
fData->Reset();
THnSparse *h1 = aGrid1->GetGrid();
THnSparse *h2 = aGrid2->GetGrid();
h2->Multiply(h1);
h2->Scale(c1*c2);
fData->Add(h2);
}
void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid, Double_t c)
{
if (aGrid->GetNVar() != GetNVar()) {
AliError("Different number of variables, cannot divide the grids");
return;
}
if (!fSumW2 && aGrid->GetSumW2()) SumW2();
THnSparse *h1 = aGrid->GetGrid();
THnSparse *h2 = (THnSparse*)fData->Clone();
fData->Divide(h2,h1);
fData->Scale(c);
}
void AliCFGridSparse::Divide(const AliCFGridSparse* aGrid1, const AliCFGridSparse* aGrid2, Double_t c1,Double_t c2, Option_t *option)
{
if (GetNVar() != aGrid1->GetNVar() || GetNVar() != aGrid2->GetNVar()) {
AliError("Different number of variables, cannot divide the grids");
return;
}
if (!fSumW2 && (aGrid1->GetSumW2() || aGrid2->GetSumW2())) SumW2();
THnSparse *h1= aGrid1->GetGrid();
THnSparse *h2= aGrid2->GetGrid();
fData->Divide(h1,h2,c1,c2,option);
}
void AliCFGridSparse::Rebin(const Int_t* group)
{
for(Int_t i=0;i<GetNVar();i++){
if (group[i]!=1) AliInfo(Form(" merging bins along dimension %i in groups of %i bins", i,group[i]));
}
THnSparse *rebinned =fData->Rebin(group);
fData->Reset();
fData = rebinned;
}
void AliCFGridSparse::Scale(Long_t index, const Double_t *fact)
{
if (GetElement(index)==0 || fact[0]==0) return;
Double_t in[2], out[2];
in[0]=GetElement(index);
in[1]=GetElementError(index);
GetScaledValues(fact,in,out);
SetElement(index,out[0]);
if (fSumW2) SetElementError(index,out[1]);
}
void AliCFGridSparse::Scale(const Int_t *bin, const Double_t *fact)
{
if(GetElement(bin)==0 || fact[0]==0)return;
Double_t in[2], out[2];
in[0]=GetElement(bin);
in[1]=GetElementError(bin);
GetScaledValues(fact,in,out);
SetElement(bin,out[0]);
if(fSumW2)SetElementError(bin,out[1]);
}
void AliCFGridSparse::Scale(const Double_t *var, const Double_t *fact)
{
if(GetElement(var)==0 || fact[0]==0)return;
Double_t in[2], out[2];
in[0]=GetElement(var);
in[1]=GetElementError(var);
GetScaledValues(fact,in,out);
SetElement(var,out[0]);
if(fSumW2)SetElementError(var,out[1]);
}
void AliCFGridSparse::Scale(const Double_t* fact)
{
for (Long_t iel=0; iel<GetNFilledBins(); iel++) {
Scale(iel,fact);
}
}
Long_t AliCFGridSparse::GetEmptyBins() const {
return (GetNBinsTotal() - GetNFilledBins()) ;
}
Int_t AliCFGridSparse::CheckStats(Double_t thr) const
{
Int_t ncellsLow=0;
for (Int_t i=0; i<GetNBinsTotal(); i++) {
if (GetElement(i)<thr) ncellsLow++;
}
return ncellsLow;
}
Double_t AliCFGridSparse::GetIntegral() const
{
return fData->ComputeIntegral();
}
Long64_t AliCFGridSparse::Merge(TCollection* list)
{
if (!list)
return 0;
if (list->IsEmpty())
return 1;
TIterator* iter = list->MakeIterator();
TObject* obj;
Int_t count = 0;
while ((obj = iter->Next())) {
AliCFGridSparse* entry = dynamic_cast<AliCFGridSparse*> (obj);
if (entry == 0)
continue;
this->Add(entry);
count++;
}
return count+1;
}
void AliCFGridSparse::GetScaledValues(const Double_t *fact, const Double_t *in, Double_t *out) const{
out[0]=in[0]*fact[0];
out[1]=TMath::Sqrt(in[1]*in[1]/in[0]/in[0]
+fact[1]*fact[1]/fact[0]/fact[0])*out[0];
}
void AliCFGridSparse::Copy(TObject& c) const
{
AliCFFrame::Copy(c);
AliCFGridSparse& target = (AliCFGridSparse &) c;
target.fSumW2 = fSumW2 ;
if (fData) {
target.fData = (THnSparse*)fData->Clone();
}
}
TH1* AliCFGridSparse::Slice(Int_t iVar1, Int_t iVar2, Int_t iVar3, const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const
{
THnSparse* clone = (THnSparse*)fData->Clone();
if (varMin != 0x0 && varMax != 0x0) {
for (Int_t iAxis=0; iAxis<GetNVar(); iAxis++) SetAxisRange(clone->GetAxis(iAxis),varMin[iAxis],varMax[iAxis],useBins);
}
TH1* projection = 0x0 ;
TString name,title;
GetProjectionName (name ,iVar1,iVar2,iVar3);
GetProjectionTitle(title,iVar1,iVar2,iVar3);
if (iVar3<0) {
if (iVar2<0) {
if (iVar1 >= GetNVar() || iVar1 < 0 ) {
AliError("Non-existent variable, return NULL");
return 0x0;
}
projection = (TH1D*)clone->Projection(iVar1);
projection->SetTitle(Form("%s_proj-%s",GetTitle(),GetVarTitle(iVar1)));
for (Int_t iBin=1; iBin<=projection->GetNbinsX(); iBin++) {
Int_t origBin = GetAxis(iVar1)->GetFirst()+iBin-1;
TString binLabel = GetAxis(iVar1)->GetBinLabel(origBin) ;
if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
}
}
else {
if (iVar1 >= GetNVar() || iVar1 < 0 ||
iVar2 >= GetNVar() || iVar2 < 0 ) {
AliError("Non-existent variable, return NULL");
return 0x0;
}
projection = (TH2D*)clone->Projection(iVar2,iVar1);
for (Int_t iBin=1; iBin<=projection->GetNbinsX(); iBin++) {
Int_t origBin = GetAxis(iVar1)->GetFirst()+iBin-1;
TString binLabel = GetAxis(iVar1)->GetBinLabel(origBin) ;
if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
}
for (Int_t iBin=1; iBin<=projection->GetNbinsY(); iBin++) {
Int_t origBin = GetAxis(iVar2)->GetFirst()+iBin-1;
TString binLabel = GetAxis(iVar2)->GetBinLabel(origBin) ;
if (binLabel.CompareTo("") != 0) projection->GetYaxis()->SetBinLabel(iBin,binLabel);
}
}
}
else {
if (iVar1 >= GetNVar() || iVar1 < 0 ||
iVar2 >= GetNVar() || iVar2 < 0 ||
iVar3 >= GetNVar() || iVar3 < 0 ) {
AliError("Non-existent variable, return NULL");
return 0x0;
}
projection = (TH3D*)clone->Projection(iVar1,iVar2,iVar3);
for (Int_t iBin=1; iBin<=projection->GetNbinsX(); iBin++) {
Int_t origBin = GetAxis(iVar1)->GetFirst()+iBin-1;
TString binLabel = GetAxis(iVar1)->GetBinLabel(origBin) ;
if (binLabel.CompareTo("") != 0) projection->GetXaxis()->SetBinLabel(iBin,binLabel);
}
for (Int_t iBin=1; iBin<=projection->GetNbinsY(); iBin++) {
Int_t origBin = GetAxis(iVar2)->GetFirst()+iBin-1;
TString binLabel = GetAxis(iVar2)->GetBinLabel(origBin) ;
if (binLabel.CompareTo("") != 0) projection->GetYaxis()->SetBinLabel(iBin,binLabel);
}
for (Int_t iBin=1; iBin<=projection->GetNbinsZ(); iBin++) {
Int_t origBin = GetAxis(iVar3)->GetFirst()+iBin-1;
TString binLabel = GetAxis(iVar3)->GetBinLabel(origBin) ;
if (binLabel.CompareTo("") != 0) projection->GetZaxis()->SetBinLabel(iBin,binLabel);
}
}
projection->SetName (name .Data());
projection->SetTitle(title.Data());
delete clone;
return projection ;
}
void AliCFGridSparse::SetAxisRange(TAxis* axis, Double_t min, Double_t max, Bool_t useBins) const {
if (useBins) axis->SetRange ((Int_t)min,(Int_t)max);
else axis->SetRangeUser( min, max);
}
void AliCFGridSparse::SetRangeUser(Int_t iVar, Double_t varMin, Double_t varMax, Bool_t useBins) const {
SetAxisRange(fData->GetAxis(iVar),varMin,varMax,useBins);
TAxis* currAxis = fData->GetAxis(iVar);
TString outString = Form("%s new range: %.1f < %s < %.1f", GetName(), currAxis->GetBinLowEdge(currAxis->GetFirst()), currAxis->GetTitle(), currAxis->GetBinUpEdge(currAxis->GetLast()));
TString binLabel = currAxis->GetBinLabel(currAxis->GetFirst());
if ( ! binLabel.IsNull() ) {
outString += " ( ";
for ( Int_t ibin = currAxis->GetFirst(); ibin <= currAxis->GetLast(); ibin++ ) {
outString += Form("%s ", currAxis->GetBinLabel(ibin));
}
outString += ")";
}
AliWarning(outString.Data());
}
void AliCFGridSparse::SetRangeUser(const Double_t *varMin, const Double_t *varMax, Bool_t useBins) const {
for (Int_t iAxis=0; iAxis<GetNVar() ; iAxis++) {
SetRangeUser(iAxis,varMin[iAxis],varMax[iAxis], useBins);
}
AliInfo("AliCFGridSparse axes ranges have been modified");
}
Float_t AliCFGridSparse::GetOverFlows(Int_t ivar, Bool_t exclusive) const
{
Int_t* bin = new Int_t[GetNVar()];
memset(bin, 0, sizeof(Int_t) * GetNVar());
Float_t ovfl=0.;
for (Long64_t i = 0; i < fData->GetNbins(); i++) {
Double_t v = fData->GetBinContent(i, bin);
Bool_t add=kTRUE;
if (exclusive) {
for(Int_t j=0;j<GetNVar();j++){
if(ivar==j)continue;
if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
}
}
if(bin[ivar]==GetNBins(ivar)+1 && add) ovfl+=v;
}
delete[] bin;
return ovfl;
}
Float_t AliCFGridSparse::GetUnderFlows(Int_t ivar, Bool_t exclusive) const
{
Int_t* bin = new Int_t[GetNVar()];
memset(bin, 0, sizeof(Int_t) * GetNVar());
Float_t unfl=0.;
for (Long64_t i = 0; i < fData->GetNbins(); i++) {
Double_t v = fData->GetBinContent(i, bin);
Bool_t add=kTRUE;
if (exclusive) {
for(Int_t j=0;j<GetNVar();j++){
if(ivar==j)continue;
if((bin[j]==0) || (bin[j]==GetNBins(j)+1))add=kFALSE;
}
}
if(bin[ivar]==0 && add) unfl+=v;
}
delete[] bin;
return unfl;
}
void AliCFGridSparse::Smooth() {
AliInfo("Your GridSparse is going to be smoothed");
AliInfo(Form("N TOTAL BINS : %li",GetNBinsTotal()));
AliInfo(Form("N FILLED BINS : %li",GetNFilledBins()));
AliCFUnfolding::SmoothUsingNeighbours(fData);
}