[u/mrichter/AliRoot.git] / PWG / Tools / AliTHn.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: AliTHn.cxx 20164 2007-08-14 15:31:50Z morsch $ */

// Use AliTHn instead of AliCFContainer and your memory consumption will be drastically reduced
// As AliTHn derives from AliCFContainer, you can just replace your current AliCFContainer object by AliTHn
// Once you have the merged output, call FillParent() and you can use AliCFContainer as usual
//
// this storage container is optimized for small memory usage
//   under/over flow bins do not exist
//   sumw2 structure is float only and only create when the weight != 1
// all histogram functionality (projections, axis, ranges, etc) are taken from THnSparse by propagating 
// the information up into the THnSparse structure (in the ananalysis *after* data processing and merging)
//
// the derivation from THnSparse is obviously against many OO rules. correct would be a common baseclass of THnSparse and THn.
// 
// Author: Jan Fiete Grosse-Oetringhaus

#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)
{
  // Constructor
}

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)
{
  // Constructor

  fNBins = 1;
  for (Int_t i=0; i<fNVars; i++)
    fNBins *= nBinIn[i];
  
  Init();
}

void AliTHn::Init()
{
  // initialize
  
  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)
{
  //
  // AliTHn copy constructor
  //

  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()
{
  // Destructor
  
  DeleteContainers();
  
  delete[] fValues;
  delete[] fSumw2;
  delete[] axisCache;
  delete[] fNbinsCache;
  delete[] fLastVars;
  delete[] fLastBins;
}

void AliTHn::DeleteContainers()
{
  // delete data containers
  
  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)
{
  // assigment operator

  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
{
  // copy function

  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)
{
  // Merge a list of AliTHn objects with this (needed for
  // PROOF). 
  // Returns the number of merged objects (including this).

  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)
{
  // fills an entry

  // fill axis cache
  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];
    
    // initial values to prevent checking for 0 below
    for (Int_t i=0; i<fNVars; i++)
    {
      fLastBins[i] = axisCache[i]->FindBin(var[i]);
      fLastVars[i] = var[i];
    }
  }
  
  // calculate global bin index
  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];
    }
    //Printf("%d", tmpBin);

    // under/overflow not supported
    if (tmpBin < 1 || tmpBin > fNbinsCache[i])
      return;
    
    // bins start from 0 here
    bin += tmpBin - 1;
//     Printf("%lld", bin);
  }

  if (!fValues[istep])
  {
    fValues[istep] = new TArrayF(fNBins);
    AliInfo(Form("Created values container for step %d", istep));
  }

  if (weight != 1)
  {
    // initialize with already filled entries (which have been filled with weight == 1), in this case fSumw2 := fValues
    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;
  
//   Printf("%f", fValues[istep][bin]);
  
  // debug
//   AliCFContainer::Fill(var, istep, weight);
}

Long64_t AliTHn::GetGlobalBinIndex(const Int_t* binIdx)
{
  // calculates global bin index
  // binIdx contains TAxis bin indexes
  // here bin count starts at 0 because we do not have over/underflow bins
  
  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)
{
  // fills the information stored in the buffer in this class into the container <cont>
  
  for (Int_t i=0; i<fNSteps; i++)
  {
    if (!fValues[i])
      continue;
      
    Float_t* source = fValues[i]->GetArray();
    // if fSumw2 is not stored, the sqrt of the number of bin entries in source is filled below; otherwise we use fSumw2
    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)
    {
//       for (Int_t j=0; j<fNVars; j++)
// 	printf("%d ", binIdx[j]);
      
      Long64_t globalBin = GetGlobalBinIndex(binIdx);
//       Printf(" --> %lld", globalBin);
      
      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()
{
  // fills the information stored in the buffer in this class into the baseclass containers
  
  FillContainer(this);
}

void AliTHn::ReduceAxis()
{
  // "removes" one axis by summing over the axis and putting the entry to bin 1
  // TODO presently only implemented for the last axis
  
  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)
    {
      // sum over axis <axis>
      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++;

      // next bin
      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;
  }
}
Commit	Line	Data
85bfac17	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	/* $Id: AliTHn.cxx 20164 2007-08-14 15:31:50Z morsch $ */
	17
f651ecbe	18	// Use AliTHn instead of AliCFContainer and your memory consumption will be drastically reduced
	19	// As AliTHn derives from AliCFContainer, you can just replace your current AliCFContainer object by AliTHn
	20	// Once you have the merged output, call FillParent() and you can use AliCFContainer as usual
85bfac17	21	//
	22	// this storage container is optimized for small memory usage
	23	// under/over flow bins do not exist
	24	// sumw2 structure is float only and only create when the weight != 1
	25	// all histogram functionality (projections, axis, ranges, etc) are taken from THnSparse by propagating
	26	// the information up into the THnSparse structure (in the ananalysis after data processing and merging)
	27	//
	28	// the derivation from THnSparse is obviously against many OO rules. correct would be a common baseclass of THnSparse and THn.
	29	//
	30	// Author: Jan Fiete Grosse-Oetringhaus
	31
	32	#include "AliTHn.h"
	33	#include "TList.h"
	34	#include "TCollection.h"
	35	#include "AliLog.h"
	36	#include "TArrayF.h"
	37	#include "THnSparse.h"
	38	#include "TMath.h"
	39
	40	ClassImp(AliTHn)
	41
	42	AliTHn::AliTHn() :
	43	AliCFContainer(),
	44	fNBins(0),
	45	fNVars(0),
	46	fNSteps(0),
	47	fValues(0),
eed401dc	48	fSumw2(0),
e97a861c	49	axisCache(0),
	50	fNbinsCache(0),
	51	fLastVars(0),
	52	fLastBins(0)
85bfac17	53	{
	54	// Constructor
	55	}
	56
	57	AliTHn::AliTHn(const Char_t* name, const Char_t* title,const Int_t nSelStep, const Int_t nVarIn, const Int_t* nBinIn) :
	58	AliCFContainer(name, title, nSelStep, nVarIn, nBinIn),
	59	fNBins(0),
	60	fNVars(nVarIn),
	61	fNSteps(nSelStep),
	62	fValues(0),
eed401dc	63	fSumw2(0),
e97a861c	64	axisCache(0),
	65	fNbinsCache(0),
	66	fLastVars(0),
	67	fLastBins(0)
85bfac17	68	{
	69	// Constructor
	70
	71	fNBins = 1;
	72	for (Int_t i=0; i<fNVars; i++)
	73	fNBins *= nBinIn[i];
	74
	75	Init();
	76	}
	77
	78	void AliTHn::Init()
	79	{
	80	// initialize
	81
	82	fValues = new TArrayF*[fNSteps];
	83	fSumw2 = new TArrayF*[fNSteps];
	84
	85	for (Int_t i=0; i<fNSteps; i++)
	86	{
	87	fValues[i] = 0;
	88	fSumw2[i] = 0;
	89	}
	90	}
	91
	92	AliTHn::AliTHn(const AliTHn &c) :
03631860	93	AliCFContainer(c),
	94	fNBins(c.fNBins),
	95	fNVars(c.fNVars),
	96	fNSteps(c.fNSteps),
	97	fValues(new TArrayF*[c.fNSteps]),
	98	fSumw2(new TArrayF*[c.fNSteps]),
e97a861c	99	axisCache(0),
	100	fNbinsCache(0),
	101	fLastVars(0),
	102	fLastBins(0)
85bfac17	103	{
	104	//
	105	// AliTHn copy constructor
	106	//
	107
03631860	108	memset(fValues,0,fNStepssizeof(TArrayF));
	109	memset(fSumw2,0,fNStepssizeof(TArrayF));
	110
	111	for (Int_t i=0; i<fNSteps; i++) {
	112	if (c.fValues[i]) fValues[i] = new TArrayF(*(c.fValues[i]));
	113	if (c.fSumw2[i]) fSumw2[i] = new TArrayF(*(c.fSumw2[i]));
	114	}
	115
85bfac17	116	}
	117
	118	AliTHn::~AliTHn()
	119	{
	120	// Destructor
	121
	122	DeleteContainers();
	123
03631860	124	delete[] fValues;
	125	delete[] fSumw2;
	126	delete[] axisCache;
e97a861c	127	delete[] fNbinsCache;
	128	delete[] fLastVars;
	129	delete[] fLastBins;
85bfac17	130	}
	131
	132	void AliTHn::DeleteContainers()
	133	{
	134	// delete data containers
	135
	136	for (Int_t i=0; i<fNSteps; i++)
	137	{
	138	if (fValues && fValues[i])
	139	{
	140	delete fValues[i];
	141	fValues[i] = 0;
	142	}
	143
	144	if (fSumw2 && fSumw2[i])
	145	{
	146	delete fSumw2[i];
	147	fSumw2[i] = 0;
	148	}
	149	}
	150	}
	151
	152	//____________________________________________________________________
	153	AliTHn &AliTHn::operator=(const AliTHn &c)
	154	{
	155	// assigment operator
	156
03631860	157	if (this != &c) {
	158	AliCFContainer::operator=(c);
	159	fNBins=c.fNBins;
	160	fNVars=c.fNVars;
03631860	161	if(fNSteps) {
	162	for(Int_t i=0; i< fNSteps; ++i) {
	163	delete fValues[i];
	164	delete fSumw2[i];
	165	}
	166	delete [] fValues;
	167	delete [] fSumw2;
3e1096ac	168	}
	169	fNSteps=c.fNSteps;
	170	if(fNSteps) {
03631860	171	fValues=new TArrayF*[fNSteps];
	172	fSumw2=new TArrayF*[fNSteps];
	173	memset(fValues,0,fNStepssizeof(TArrayF));
	174	memset(fSumw2,0,fNStepssizeof(TArrayF));
	175
	176	for (Int_t i=0; i<fNSteps; i++) {
	177	if (c.fValues[i]) fValues[i] = new TArrayF(*(c.fValues[i]));
	178	if (c.fSumw2[i]) fSumw2[i] = new TArrayF(*(c.fSumw2[i]));
	179	}
3e1096ac	180	} else {
	181	fValues = 0;
	182	fSumw2 = 0;
03631860	183	}
03631860	184	delete [] axisCache;
32aa0326	185	axisCache = new TAxis*[fNVars];
25a9e954	186	memcpy(axisCache, c.axisCache, fNVarssizeof(TAxis));
03631860	187	}
85bfac17	188	return *this;
	189	}
	190
	191	//____________________________________________________________________
	192	void AliTHn::Copy(TObject& c) const
	193	{
	194	// copy function
	195
	196	AliTHn& target = (AliTHn &) c;
	197
	198	AliCFContainer::Copy(target);
	199
	200	target.fNSteps = fNSteps;
	201	target.fNBins = fNBins;
	202	target.fNVars = fNVars;
	203
	204	target.Init();
	205
	206	for (Int_t i=0; i<fNSteps; i++)
	207	{
	208	if (fValues[i])
	209	target.fValues[i] = new TArrayF(*(fValues[i]));
	210	else
	211	target.fValues[i] = 0;
	212
	213	if (fSumw2[i])
	214	target.fSumw2[i] = new TArrayF(*(fSumw2[i]));
	215	else
	216	target.fSumw2[i] = 0;
	217	}
	218	}
	219
	220	//____________________________________________________________________
	221	Long64_t AliTHn::Merge(TCollection* list)
	222	{
	223	// Merge a list of AliTHn objects with this (needed for
	224	// PROOF).
	225	// Returns the number of merged objects (including this).
	226
	227	if (!list)
	228	return 0;
	229
	230	if (list->IsEmpty())
	231	return 1;
	232
	233	AliCFContainer::Merge(list);
	234
	235	TIterator* iter = list->MakeIterator();
	236	TObject* obj;
	237
	238	Int_t count = 0;
	239	while ((obj = iter->Next())) {
	240
	241	AliTHn* entry = dynamic_cast<AliTHn*> (obj);
	242	if (entry == 0)
	243	continue;
	244
	245	for (Int_t i=0; i<fNSteps; i++)
	246	{
	247	if (entry->fValues[i])
	248	{
	249	if (!fValues[i])
	250	fValues[i] = new TArrayF(fNBins);
	251
252	for (Long64_t l = 0; l<fNBins; l++)
253	fValues[i]->GetArray()[l] += entry->fValues[i]->GetArray()[l];
254	}
255
256	if (entry->fSumw2[i])
257	{
258	if (!fSumw2[i])
259	fSumw2[i] = new TArrayF(fNBins);
260
261	for (Long64_t l = 0; l<fNBins; l++)
262	fSumw2[i]->GetArray()[l] += entry->fSumw2[i]->GetArray()[l];
263	}
264	}
265
266	count++;
267	}
268
269	return count+1;
270	}
271
272	void AliTHn::Fill(const Double_t *var, Int_t istep, Double_t weight)
273	{
274	// fills an entry
275
eed401dc	276	// fill axis cache
	277	if (!axisCache)
	278	{
	279	axisCache = new TAxis*[fNVars];
e97a861c	280	fNbinsCache = new Int_t[fNVars];
eed401dc	281	for (Int_t i=0; i<fNVars; i++)
e97a861c	282	{
eed401dc	283	axisCache[i] = GetAxis(i, 0);
e97a861c	284	fNbinsCache[i] = axisCache[i]->GetNbins();
	285	}
	286
	287	fLastVars = new Double_t[fNVars];
	288	fLastBins = new Int_t[fNVars];
	289
	290	// initial values to prevent checking for 0 below
	291	for (Int_t i=0; i<fNVars; i++)
	292	{
	293	fLastBins[i] = axisCache[i]->FindBin(var[i]);
	294	fLastVars[i] = var[i];
	295	}
eed401dc	296	}
e97a861c	297
85bfac17	298	// calculate global bin index
	299	Long64_t bin = 0;
	300	for (Int_t i=0; i<fNVars; i++)
	301	{
e97a861c	302	bin *= fNbinsCache[i];
85bfac17	303
e97a861c	304	Int_t tmpBin = 0;
	305	if (fLastVars[i] == var[i])
	306	tmpBin = fLastBins[i];
	307	else
	308	{
	309	tmpBin = axisCache[i]->FindBin(var[i]);
	310	fLastBins[i] = tmpBin;
	311	fLastVars[i] = var[i];
	312	}
	313	//Printf("%d", tmpBin);
	314
85bfac17	315	// under/overflow not supported
e97a861c	316	if (tmpBin < 1 \|\| tmpBin > fNbinsCache[i])
85bfac17	317	return;
	318
	319	// bins start from 0 here
	320	bin += tmpBin - 1;
	321	// Printf("%lld", bin);
	322	}
	323
	324	if (!fValues[istep])
	325	{
	326	fValues[istep] = new TArrayF(fNBins);
	327	AliInfo(Form("Created values container for step %d", istep));
	328	}
	329
	330	if (weight != 1)
	331	{
	332	// initialize with already filled entries (which have been filled with weight == 1), in this case fSumw2 := fValues
	333	if (!fSumw2[istep])
	334	{
	335	fSumw2[istep] = new TArrayF(*fValues[istep]);
	336	AliInfo(Form("Created sumw2 container for step %d", istep));
	337	}
	338	}
	339
	340	fValues[istep]->GetArray()[bin] += weight;
	341	if (fSumw2[istep])
	342	fSumw2[istep]->GetArray()[bin] += weight * weight;
	343
	344	// Printf("%f", fValues[istep][bin]);
	345
	346	// debug
	347	// AliCFContainer::Fill(var, istep, weight);
	348	}
	349
	350	Long64_t AliTHn::GetGlobalBinIndex(const Int_t* binIdx)
	351	{
	352	// calculates global bin index
	353	// binIdx contains TAxis bin indexes
	354	// here bin count starts at 0 because we do not have over/underflow bins
	355
	356	Long64_t bin = 0;
	357	for (Int_t i=0; i<fNVars; i++)
	358	{
	359	bin *= GetAxis(i, 0)->GetNbins();
	360	bin += binIdx[i] - 1;
	361	}
	362
	363	return bin;
	364	}
	365
1bba939a	366	void AliTHn::FillContainer(AliCFContainer* cont)
85bfac17	367	{
1bba939a	368	// fills the information stored in the buffer in this class into the container <cont>
85bfac17	369
	370	for (Int_t i=0; i<fNSteps; i++)
	371	{
	372	if (!fValues[i])
	373	continue;
	374
	375	Float_t* source = fValues[i]->GetArray();
	376	// if fSumw2 is not stored, the sqrt of the number of bin entries in source is filled below; otherwise we use fSumw2
	377	Float_t* sourceSumw2 = source;
	378	if (fSumw2[i])
	379	sourceSumw2 = fSumw2[i]->GetArray();
	380
1bba939a	381	THnSparse* target = cont->GetGrid(i)->GetGrid();
85bfac17	382
85bfac17	383	Int_t* binIdx = new Int_t[fNVars];
c32a0ca9	384	Int_t* nBins = new Int_t[fNVars];
85bfac17	385	for (Int_t j=0; j<fNVars; j++)
c32a0ca9	386	{
85bfac17	387	binIdx[j] = 1;
c32a0ca9	388	nBins[j] = target->GetAxis(j)->GetNbins();
c32a0ca9	389	}
85bfac17	390
	391	Long64_t count = 0;
	392
	393	while (1)
	394	{
	395	// for (Int_t j=0; j<fNVars; j++)
	396	// printf("%d ", binIdx[j]);
	397
	398	Long64_t globalBin = GetGlobalBinIndex(binIdx);
	399	// Printf(" --> %lld", globalBin);
	400
	401	if (source[globalBin] != 0)
	402	{
	403	target->SetBinContent(binIdx, source[globalBin]);
	404	target->SetBinError(binIdx, TMath::Sqrt(sourceSumw2[globalBin]));
	405
	406	count++;
	407	}
	408
	409	binIdx[fNVars-1]++;
	410
	411	for (Int_t j=fNVars-1; j>0; j--)
	412	{
c32a0ca9	413	if (binIdx[j] > nBins[j])
85bfac17	414	{
	415	binIdx[j] = 1;
	416	binIdx[j-1]++;
	417	}
	418	}
	419
c32a0ca9	420	if (binIdx[0] > nBins[0])
85bfac17	421	break;
	422	}
	423
	424	AliInfo(Form("Step %d: copied %lld entries out of %lld bins", i, count, GetGlobalBinIndex(binIdx)));
c32a0ca9	425
	426	delete[] binIdx;
	427	delete[] nBins;
1bba939a	428	}
	429	}
	430
	431	void AliTHn::FillParent()
	432	{
	433	// fills the information stored in the buffer in this class into the baseclass containers
	434
	435	FillContainer(this);
c32a0ca9	436	}
	437
	438	void AliTHn::ReduceAxis()
	439	{
	440	// "removes" one axis by summing over the axis and putting the entry to bin 1
	441	// TODO presently only implemented for the last axis
	442
	443	Int_t axis = fNVars-1;
	444
	445	for (Int_t i=0; i<fNSteps; i++)
	446	{
	447	if (!fValues[i])
	448	continue;
	449
	450	Float_t* source = fValues[i]->GetArray();
	451	Float_t* sourceSumw2 = 0;
	452	if (fSumw2[i])
	453	sourceSumw2 = fSumw2[i]->GetArray();
	454
	455	THnSparse* target = GetGrid(i)->GetGrid();
	456
	457	Int_t* binIdx = new Int_t[fNVars];
	458	Int_t* nBins = new Int_t[fNVars];
	459	for (Int_t j=0; j<fNVars; j++)
	460	{
	461	binIdx[j] = 1;
	462	nBins[j] = target->GetAxis(j)->GetNbins();
	463	}
	464
	465	Long64_t count = 0;
	466
	467	while (1)
	468	{
	469	// sum over axis <axis>
	470	Float_t sumValues = 0;
	471	Float_t sumSumw2 = 0;
	472	for (Int_t j=1; j<=nBins[axis]; j++)
	473	{
	474	binIdx[axis] = j;
	475	Long64_t globalBin = GetGlobalBinIndex(binIdx);
	476	sumValues += source[globalBin];
	477	source[globalBin] = 0;
	478
	479	if (sourceSumw2)
	480	{
	481	sumSumw2 += sourceSumw2[globalBin];
	482	sourceSumw2[globalBin] = 0;
	483	}
	484	}
	485	binIdx[axis] = 1;
	486
	487	Long64_t globalBin = GetGlobalBinIndex(binIdx);
	488	source[globalBin] = sumValues;
	489	if (sourceSumw2)
	490	sourceSumw2[globalBin] = sumSumw2;
	491
	492	count++;
	493
	494	// next bin
	495	binIdx[fNVars-2]++;
	496
	497	for (Int_t j=fNVars-2; j>0; j--)
	498	{
	499	if (binIdx[j] > nBins[j])
500	{
501	binIdx[j] = 1;
502	binIdx[j-1]++;
503	}
504	}
505
506	if (binIdx[0] > nBins[0])
507	break;
508	}
509
510	AliInfo(Form("Step %d: reduced %lld bins to %lld entries", i, GetGlobalBinIndex(binIdx), count));
511
512	delete[] binIdx;
513	delete[] nBins;
85bfac17	514	}
85bfac17	515	}