[u/mrichter/AliRoot.git] / PWG0 / AliCorrectionMatrix3D.cxx

/* $Id$ */

// ------------------------------------------------------
//
// Class to handle 3d-corrections.
//
// ------------------------------------------------------
//

#include <TDirectory.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TString.h>
#include <TMath.h>

#include <AliLog.h>

#include "AliCorrectionMatrix3D.h"
#include "AliCorrectionMatrix2D.h"
#include "AliPWG0Helper.h"

//____________________________________________________________________
ClassImp(AliCorrectionMatrix3D)

//____________________________________________________________________
AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
  AliCorrectionMatrix()
{
  // default constructor
}

//____________________________________________________________________
AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
  : AliCorrectionMatrix(c)
{
  // copy constructor
  ((AliCorrectionMatrix3D &)c).Copy(*this);
}

//____________________________________________________________________
AliCorrectionMatrix3D &AliCorrectionMatrix3D::operator=(const AliCorrectionMatrix3D &c)
{
  // assigment operator

  if (this != &c)
    ((AliCorrectionMatrix3D &) c).Copy(*this);

  return *this;
}

//____________________________________________________________________
AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
              Int_t nBinX, Float_t Xmin, Float_t Xmax,
              Int_t nBinY, Float_t Ymin, Float_t Ymax,
              Int_t nBinZ, Float_t Zmin, Float_t Zmax)
  : AliCorrectionMatrix(name, title)
{
  //
  // constructor
  //

  Float_t* binLimitsX = new Float_t[nBinX+1];
  for (Int_t i=0; i<=nBinX; ++i)
    binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;

  Float_t* binLimitsY = new Float_t[nBinY+1];
  for (Int_t i=0; i<=nBinY; ++i)
    binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;

  Float_t* binLimitsZ = new Float_t[nBinZ+1];
  for (Int_t i=0; i<=nBinZ; ++i)
    binLimitsZ[i] = Zmin + (Zmax - Zmin) / nBinZ * i;

  CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);

  delete[] binLimitsX;
  delete[] binLimitsY;
  delete[] binLimitsZ;
}

AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
      Int_t nBinX, Float_t Xmin, Float_t Xmax,
      Int_t nBinY, Float_t Ymin, Float_t Ymax,
      Int_t nBinZ, const Float_t* zbins)
  : AliCorrectionMatrix(name, title)
{
  // constructor with variable bin sizes

  Float_t* binLimitsX = new Float_t[nBinX+1];
  for (Int_t i=0; i<=nBinX; ++i)
    binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;

  Float_t* binLimitsY = new Float_t[nBinY+1];
  for (Int_t i=0; i<=nBinY; ++i)
    binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;

  CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);

  delete[] binLimitsX;
  delete[] binLimitsY;
}

AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3* hBinning)
  : AliCorrectionMatrix(name, title)
{
  // constructor with variable bin sizes (uses binning of hBinning)

  // do not add this hists to the directory
  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  fhMeas = (TH3F*)hBinning->Clone("measured");
  fhGene = (TH3F*)hBinning->Clone("generated");
  fhCorr = (TH3F*)hBinning->Clone("correction");

  fhMeas->SetTitle(Form("%s measured", GetTitle()));
  fhGene->SetTitle(Form("%s generated", GetTitle()));
  fhCorr->SetTitle(Form("%s correction", GetTitle()));

  fhMeas->Reset();
  fhGene->Reset();
  fhCorr->Reset();

  TH1::AddDirectory(oldStatus);

  fhMeas->Sumw2();
  fhGene->Sumw2();
  fhCorr->Sumw2();
}

//____________________________________________________________________
void AliCorrectionMatrix3D::CreateHists(Int_t nBinX, const Float_t* binLimitsX,
      Int_t nBinY, const Float_t* binLimitsY,
      Int_t nBinZ, const Float_t* binLimitsZ)
{
  // create the histograms

  // do not add this hists to the directory
  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  fhMeas  = new TH3F("measured",   Form("%s measured",GetTitle()),   nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
  fhGene  = new TH3F("generated",  Form("%s generated",GetTitle()),  nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
  fhCorr  = new TH3F("correction", Form("%s correction",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);

  fhMeas->Sumw2();
  fhGene->Sumw2();
  fhCorr->Sumw2();

  TH1::AddDirectory(oldStatus);
}


//____________________________________________________________________
AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
{
  //
  // destructor
  //

  // histograms already deleted in base class
}

//____________________________________________________________________
TH3* AliCorrectionMatrix3D::GetGeneratedHistogram()
{
  // return generated histogram casted to correct type
  return dynamic_cast<TH3F*> (fhGene);
}

//____________________________________________________________________
TH3* AliCorrectionMatrix3D::GetMeasuredHistogram()
{
  // return measured histogram casted to correct type
  return dynamic_cast<TH3F*> (fhMeas);
}

//____________________________________________________________________
TH3* AliCorrectionMatrix3D::GetCorrectionHistogram()
{
  // return correction histogram casted to correct type
  return dynamic_cast<TH3F*> (fhCorr);
}

//____________________________________________________________________
AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Option_t* opt, Float_t aMin, Float_t aMax)
{
  // returns a 2D projection of this correction

  TString option = opt;

  // unzoom
  fhMeas->GetXaxis()->SetRange(0, 0);
  fhMeas->GetYaxis()->SetRange(0, 0);
  fhMeas->GetZaxis()->SetRange(0, 0);

  fhGene->GetXaxis()->SetRange(0, 0);
  fhGene->GetYaxis()->SetRange(0, 0);
  fhGene->GetZaxis()->SetRange(0, 0);

  if (aMin<aMax) {
    if (option.Contains("xy") || option.Contains("yx")) {
      Int_t bMin = fhMeas->GetZaxis()->FindBin(aMin);
      Int_t bMax = fhMeas->GetZaxis()->FindBin(aMax);
      fhGene->GetZaxis()->SetRange(bMin, bMax);
      fhMeas->GetZaxis()->SetRange(bMin, bMax);
    }
    else if (option.Contains("xz") || option.Contains("zx")) {
      Int_t bMin = fhMeas->GetYaxis()->FindBin(aMin);
      Int_t bMax = fhMeas->GetYaxis()->FindBin(aMax);
      fhGene->GetYaxis()->SetRange(bMin, bMax);
      fhMeas->GetYaxis()->SetRange(bMin, bMax);
    }
    else if (option.Contains("yz") || option.Contains("zy")) {
      Int_t bMin = fhMeas->GetXaxis()->FindBin(aMin);
      Int_t bMax = fhMeas->GetXaxis()->FindBin(aMax);
      fhGene->GetXaxis()->SetRange(bMin, bMax);
      fhMeas->GetXaxis()->SetRange(bMin, bMax);
    }
    else {
      AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
      return 0;
    }
  }

  AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);

  TH2F* meas = (TH2F*) ((TH3F*)fhMeas)->Project3D(option)->Clone(Form("%s_meas", corr2D->GetName()));
  TH2F* gene = (TH2F*) ((TH3F*)fhGene)->Project3D(option)->Clone(Form("%s_gene", corr2D->GetName()));
  
  // set errors
  for (Int_t x = 0; x<=meas->GetNbinsX()+1; x++)
    for (Int_t y = 0; y<=meas->GetNbinsY()+1; y++)
    {
      gene->SetBinError(x, y, TMath::Sqrt(gene->GetBinContent(x, y)));
      meas->SetBinError(x, y, TMath::Sqrt(meas->GetBinContent(x, y)));
    }

  TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
  corr->Reset();

  corr2D->SetGeneratedHistogram(gene);
  corr2D->SetMeasuredHistogram(meas);
  corr2D->SetCorrectionHistogram(corr);

  corr2D->Divide();

  // unzoom
  fhMeas->GetXaxis()->SetRange(0, 0);
  fhMeas->GetYaxis()->SetRange(0, 0);
  fhMeas->GetZaxis()->SetRange(0, 0);

  fhGene->GetXaxis()->SetRange(0, 0);
  fhGene->GetYaxis()->SetRange(0, 0);
  fhGene->GetZaxis()->SetRange(0, 0);

  return corr2D;
}

//____________________________________________________________________
TH1* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Option_t* opt, Float_t aMin1, Float_t aMax1, Float_t aMin2, Float_t aMax2)
{
  // returns a 1D projection of this correction
  
  AliCorrectionMatrix2D* corr2D = 0;
  if (strcmp(opt,"x")==0) {
    corr2D = Get2DCorrection("yx",aMin2,aMax2);
    return corr2D->Get1DCorrectionHistogram("x",aMin1,aMax1);
  }
  if (strcmp(opt,"y")==0) {
    corr2D = Get2DCorrection("xy",aMin2,aMax2);
    return corr2D->Get1DCorrectionHistogram("x",aMin1,aMax1);
  }  
  if (strcmp(opt,"z")==0) {
    corr2D = Get2DCorrection("yz",aMin1,aMax1);
    return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
  }  
  AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));
  
  return 0;
}

//____________________________________________________________________
void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az, Double_t weight)
{
  // add value to measured histogram
  GetMeasuredHistogram()->Fill(ax, ay, az, weight);
}

//____________________________________________________________________
void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az, Double_t weight)
{
  // add value to generated histogram
  GetGeneratedHistogram()->Fill(ax, ay, az, weight);
}

//____________________________________________________________________
Float_t AliCorrectionMatrix3D::GetCorrection(Float_t ax, Float_t ay, Float_t az) const
{
  // returns a value of the correction map
  return fhCorr->GetBinContent(fhCorr->FindBin(ax, ay, az));
}

//____________________________________________________________________
//void AliCorrectionMatrix3D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge, Int_t nBinsZedge)
void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
{
  // so what do we do here...
}

//____________________________________________________________________
void AliCorrectionMatrix3D::SaveHistograms()
{
  //
  // saves the histograms
  //

  AliCorrectionMatrix::SaveHistograms();

  if (GetGeneratedHistogram() && GetMeasuredHistogram())
  {
    gDirectory->cd(GetName());
    
    AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
    
    gDirectory->cd("..");
  }
}

//____________________________________________________________________
Int_t AliCorrectionMatrix3D::CheckEmptyBins(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax, Bool_t quiet)
{
  //
  // counts the number of empty Bins in a given region
  //

  TH3* hist = GetGeneratedHistogram();
  if (!hist)
    return -1;

  Int_t emptyBins = 0;
  for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
    for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
      for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
        if (hist->GetBinContent(x, y, z) == 0)
        {
          if (!quiet)
            printf("Empty bin in %s at vtx = %f, eta = %f, pt = %f\n", GetName(), hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
          ++emptyBins;
        }

  return emptyBins;
}

//____________________________________________________________________
TH1* AliCorrectionMatrix3D::PlotBinErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax)
{
  //
  // makes a 1d plots of the relative bin errors
  //

  TH3* hist = GetCorrectionHistogram();
  if (!hist)
    return 0;

  TH1F* target = new TH1F("relerrors", "relerrors", 100, 0, 10);

  for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
    for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
      for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
        if (hist->GetBinContent(x, y, z) != 0)
          target->Fill(100 * hist->GetBinError(x, y, z) / hist->GetBinContent(x, y, z));

  return target;
}
Commit	Line	Data
	1	/* $Id$ */
	2
	3	// ------------------------------------------------------
	4	//
	5	// Class to handle 3d-corrections.
	6	//
	7	// ------------------------------------------------------
	8	//
	9
	10	#include <TDirectory.h>
	11	#include <TH1F.h>
	12	#include <TH2F.h>
	13	#include <TH3F.h>
	14	#include <TString.h>
	15	#include <TMath.h>
	16
	17	#include <AliLog.h>
	18
	19	#include "AliCorrectionMatrix3D.h"
	20	#include "AliCorrectionMatrix2D.h"
	21	#include "AliPWG0Helper.h"
	22
	23	//____________________________________________________________________
	24	ClassImp(AliCorrectionMatrix3D)
	25
	26	//____________________________________________________________________
	27	AliCorrectionMatrix3D::AliCorrectionMatrix3D() :
	28	AliCorrectionMatrix()
	29	{
	30	// default constructor
	31	}
	32
	33	//____________________________________________________________________
	34	AliCorrectionMatrix3D::AliCorrectionMatrix3D(const AliCorrectionMatrix3D& c)
	35	: AliCorrectionMatrix(c)
	36	{
	37	// copy constructor
	38	((AliCorrectionMatrix3D &)c).Copy(*this);
	39	}
	40
	41	//____________________________________________________________________
	42	AliCorrectionMatrix3D &AliCorrectionMatrix3D::operator=(const AliCorrectionMatrix3D &c)
	43	{
	44	// assigment operator
	45
	46	if (this != &c)
	47	((AliCorrectionMatrix3D &) c).Copy(*this);
	48
	49	return *this;
	50	}
	51
	52	//____________________________________________________________________
	53	AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
	54	Int_t nBinX, Float_t Xmin, Float_t Xmax,
	55	Int_t nBinY, Float_t Ymin, Float_t Ymax,
	56	Int_t nBinZ, Float_t Zmin, Float_t Zmax)
	57	: AliCorrectionMatrix(name, title)
	58	{
	59	//
	60	// constructor
	61	//
	62
	63	Float_t* binLimitsX = new Float_t[nBinX+1];
	64	for (Int_t i=0; i<=nBinX; ++i)
	65	binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
	66
	67	Float_t* binLimitsY = new Float_t[nBinY+1];
	68	for (Int_t i=0; i<=nBinY; ++i)
	69	binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
	70
	71	Float_t* binLimitsZ = new Float_t[nBinZ+1];
	72	for (Int_t i=0; i<=nBinZ; ++i)
	73	binLimitsZ[i] = Zmin + (Zmax - Zmin) / nBinZ * i;
	74
	75	CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
	76
	77	delete[] binLimitsX;
	78	delete[] binLimitsY;
	79	delete[] binLimitsZ;
	80	}
	81
	82	AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title,
	83	Int_t nBinX, Float_t Xmin, Float_t Xmax,
	84	Int_t nBinY, Float_t Ymin, Float_t Ymax,
	85	Int_t nBinZ, const Float_t* zbins)
	86	: AliCorrectionMatrix(name, title)
	87	{
	88	// constructor with variable bin sizes
	89
	90	Float_t* binLimitsX = new Float_t[nBinX+1];
	91	for (Int_t i=0; i<=nBinX; ++i)
	92	binLimitsX[i] = Xmin + (Xmax - Xmin) / nBinX * i;
	93
	94	Float_t* binLimitsY = new Float_t[nBinY+1];
	95	for (Int_t i=0; i<=nBinY; ++i)
	96	binLimitsY[i] = Ymin + (Ymax - Ymin) / nBinY * i;
	97
	98	CreateHists(nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, zbins);
	99
	100	delete[] binLimitsX;
	101	delete[] binLimitsY;
	102	}
	103
	104	AliCorrectionMatrix3D::AliCorrectionMatrix3D(const Char_t* name, const Char_t* title, TH3* hBinning)
	105	: AliCorrectionMatrix(name, title)
	106	{
	107	// constructor with variable bin sizes (uses binning of hBinning)
	108
	109	// do not add this hists to the directory
	110	Bool_t oldStatus = TH1::AddDirectoryStatus();
	111	TH1::AddDirectory(kFALSE);
	112
	113	fhMeas = (TH3F*)hBinning->Clone("measured");
	114	fhGene = (TH3F*)hBinning->Clone("generated");
	115	fhCorr = (TH3F*)hBinning->Clone("correction");
	116
	117	fhMeas->SetTitle(Form("%s measured", GetTitle()));
	118	fhGene->SetTitle(Form("%s generated", GetTitle()));
	119	fhCorr->SetTitle(Form("%s correction", GetTitle()));
	120
	121	fhMeas->Reset();
	122	fhGene->Reset();
	123	fhCorr->Reset();
	124
	125	TH1::AddDirectory(oldStatus);
	126
	127	fhMeas->Sumw2();
	128	fhGene->Sumw2();
	129	fhCorr->Sumw2();
	130	}
	131
	132	//____________________________________________________________________
	133	void AliCorrectionMatrix3D::CreateHists(Int_t nBinX, const Float_t* binLimitsX,
	134	Int_t nBinY, const Float_t* binLimitsY,
	135	Int_t nBinZ, const Float_t* binLimitsZ)
	136	{
	137	// create the histograms
	138
	139	// do not add this hists to the directory
	140	Bool_t oldStatus = TH1::AddDirectoryStatus();
	141	TH1::AddDirectory(kFALSE);
	142
	143	fhMeas = new TH3F("measured", Form("%s measured",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
	144	fhGene = new TH3F("generated", Form("%s generated",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
	145	fhCorr = new TH3F("correction", Form("%s correction",GetTitle()), nBinX, binLimitsX, nBinY, binLimitsY, nBinZ, binLimitsZ);
	146
	147	fhMeas->Sumw2();
	148	fhGene->Sumw2();
	149	fhCorr->Sumw2();
	150
	151	TH1::AddDirectory(oldStatus);
	152	}
	153
	154
	155	//____________________________________________________________________
	156	AliCorrectionMatrix3D::~AliCorrectionMatrix3D()
	157	{
	158	//
	159	// destructor
	160	//
	161
	162	// histograms already deleted in base class
	163	}
	164
	165	//____________________________________________________________________
	166	TH3* AliCorrectionMatrix3D::GetGeneratedHistogram()
	167	{
	168	// return generated histogram casted to correct type
	169	return dynamic_cast<TH3F*> (fhGene);
	170	}
	171
	172	//____________________________________________________________________
	173	TH3* AliCorrectionMatrix3D::GetMeasuredHistogram()
	174	{
	175	// return measured histogram casted to correct type
	176	return dynamic_cast<TH3F*> (fhMeas);
	177	}
	178
	179	//____________________________________________________________________
	180	TH3* AliCorrectionMatrix3D::GetCorrectionHistogram()
	181	{
	182	// return correction histogram casted to correct type
	183	return dynamic_cast<TH3F*> (fhCorr);
	184	}
	185
	186	//____________________________________________________________________
	187	AliCorrectionMatrix2D* AliCorrectionMatrix3D::Get2DCorrection(Option_t* opt, Float_t aMin, Float_t aMax)
	188	{
	189	// returns a 2D projection of this correction
	190
	191	TString option = opt;
	192
	193	// unzoom
	194	fhMeas->GetXaxis()->SetRange(0, 0);
	195	fhMeas->GetYaxis()->SetRange(0, 0);
	196	fhMeas->GetZaxis()->SetRange(0, 0);
	197
	198	fhGene->GetXaxis()->SetRange(0, 0);
	199	fhGene->GetYaxis()->SetRange(0, 0);
	200	fhGene->GetZaxis()->SetRange(0, 0);
	201
	202	if (aMin<aMax) {
	203	if (option.Contains("xy") \|\| option.Contains("yx")) {
	204	Int_t bMin = fhMeas->GetZaxis()->FindBin(aMin);
	205	Int_t bMax = fhMeas->GetZaxis()->FindBin(aMax);
	206	fhGene->GetZaxis()->SetRange(bMin, bMax);
	207	fhMeas->GetZaxis()->SetRange(bMin, bMax);
	208	}
	209	else if (option.Contains("xz") \|\| option.Contains("zx")) {
	210	Int_t bMin = fhMeas->GetYaxis()->FindBin(aMin);
	211	Int_t bMax = fhMeas->GetYaxis()->FindBin(aMax);
	212	fhGene->GetYaxis()->SetRange(bMin, bMax);
	213	fhMeas->GetYaxis()->SetRange(bMin, bMax);
	214	}
	215	else if (option.Contains("yz") \|\| option.Contains("zy")) {
	216	Int_t bMin = fhMeas->GetXaxis()->FindBin(aMin);
	217	Int_t bMax = fhMeas->GetXaxis()->FindBin(aMax);
	218	fhGene->GetXaxis()->SetRange(bMin, bMax);
	219	fhMeas->GetXaxis()->SetRange(bMin, bMax);
	220	}
	221	else {
	222	AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s", opt));
	223	return 0;
	224	}
	225	}
	226
	227	AliCorrectionMatrix2D* corr2D = new AliCorrectionMatrix2D(Form("%s_%s",GetName(),opt),Form("%s projection %s",GetName(),opt),100,0,100,100,0,100);
	228
	229	TH2F* meas = (TH2F) ((TH3F)fhMeas)->Project3D(option)->Clone(Form("%s_meas", corr2D->GetName()));
	230	TH2F* gene = (TH2F) ((TH3F)fhGene)->Project3D(option)->Clone(Form("%s_gene", corr2D->GetName()));
	231
	232	// set errors
	233	for (Int_t x = 0; x<=meas->GetNbinsX()+1; x++)
	234	for (Int_t y = 0; y<=meas->GetNbinsY()+1; y++)
	235	{
	236	gene->SetBinError(x, y, TMath::Sqrt(gene->GetBinContent(x, y)));
	237	meas->SetBinError(x, y, TMath::Sqrt(meas->GetBinContent(x, y)));
	238	}
	239
	240	TH2F* corr = (TH2F*)gene->Clone(Form("%s_corr", corr2D->GetName()));
	241	corr->Reset();
	242
	243	corr2D->SetGeneratedHistogram(gene);
	244	corr2D->SetMeasuredHistogram(meas);
	245	corr2D->SetCorrectionHistogram(corr);
	246
	247	corr2D->Divide();
	248
	249	// unzoom
	250	fhMeas->GetXaxis()->SetRange(0, 0);
	251	fhMeas->GetYaxis()->SetRange(0, 0);
	252	fhMeas->GetZaxis()->SetRange(0, 0);
	253
	254	fhGene->GetXaxis()->SetRange(0, 0);
	255	fhGene->GetYaxis()->SetRange(0, 0);
	256	fhGene->GetZaxis()->SetRange(0, 0);
	257
	258	return corr2D;
	259	}
	260
	261	//____________________________________________________________________
	262	TH1* AliCorrectionMatrix3D::Get1DCorrectionHistogram(Option_t* opt, Float_t aMin1, Float_t aMax1, Float_t aMin2, Float_t aMax2)
	263	{
	264	// returns a 1D projection of this correction
	265
	266	AliCorrectionMatrix2D* corr2D = 0;
	267	if (strcmp(opt,"x")==0) {
	268	corr2D = Get2DCorrection("yx",aMin2,aMax2);
	269	return corr2D->Get1DCorrectionHistogram("x",aMin1,aMax1);
	270	}
	271	if (strcmp(opt,"y")==0) {
	272	corr2D = Get2DCorrection("xy",aMin2,aMax2);
	273	return corr2D->Get1DCorrectionHistogram("x",aMin1,aMax1);
	274	}
	275	if (strcmp(opt,"z")==0) {
	276	corr2D = Get2DCorrection("yz",aMin1,aMax1);
	277	return corr2D->Get1DCorrectionHistogram("x",aMin2,aMax2);
	278	}
	279	AliDebug(AliLog::kWarning, Form("WARNING: unknown projection option %s (should be x,y or z)", opt));
	280
	281	return 0;
	282	}
	283
	284	//____________________________________________________________________
	285	void AliCorrectionMatrix3D::FillMeas(Float_t ax, Float_t ay, Float_t az, Double_t weight)
	286	{
	287	// add value to measured histogram
	288	GetMeasuredHistogram()->Fill(ax, ay, az, weight);
	289	}
	290
	291	//____________________________________________________________________
	292	void AliCorrectionMatrix3D::FillGene(Float_t ax, Float_t ay, Float_t az, Double_t weight)
	293	{
	294	// add value to generated histogram
	295	GetGeneratedHistogram()->Fill(ax, ay, az, weight);
	296	}
	297
	298	//____________________________________________________________________
	299	Float_t AliCorrectionMatrix3D::GetCorrection(Float_t ax, Float_t ay, Float_t az) const
	300	{
	301	// returns a value of the correction map
	302	return fhCorr->GetBinContent(fhCorr->FindBin(ax, ay, az));
	303	}
	304
	305	//____________________________________________________________________
	306	//void AliCorrectionMatrix3D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge, Int_t nBinsZedge)
	307	void AliCorrectionMatrix3D::RemoveEdges(Float_t, Int_t, Int_t, Int_t)
	308	{
	309	// so what do we do here...
	310	}
	311
	312	//____________________________________________________________________
	313	void AliCorrectionMatrix3D::SaveHistograms()
	314	{
	315	//
	316	// saves the histograms
	317	//
	318
	319	AliCorrectionMatrix::SaveHistograms();
	320
	321	if (GetGeneratedHistogram() && GetMeasuredHistogram())
	322	{
	323	gDirectory->cd(GetName());
	324
	325	AliPWG0Helper::CreateDividedProjections(GetGeneratedHistogram(), GetMeasuredHistogram(), 0, kFALSE, kTRUE);
	326
	327	gDirectory->cd("..");
	328	}
	329	}
	330
	331	//____________________________________________________________________
	332	Int_t AliCorrectionMatrix3D::CheckEmptyBins(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax, Bool_t quiet)
	333	{
	334	//
	335	// counts the number of empty Bins in a given region
	336	//
	337
	338	TH3* hist = GetGeneratedHistogram();
	339	if (!hist)
	340	return -1;
	341
	342	Int_t emptyBins = 0;
	343	for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
	344	for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
	345	for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
	346	if (hist->GetBinContent(x, y, z) == 0)
	347	{
	348	if (!quiet)
	349	printf("Empty bin in %s at vtx = %f, eta = %f, pt = %f\n", GetName(), hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
	350	++emptyBins;
	351	}
	352
	353	return emptyBins;
	354	}
	355
	356	//____________________________________________________________________
	357	TH1* AliCorrectionMatrix3D::PlotBinErrors(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t zmin, Float_t zmax)
	358	{
	359	//
	360	// makes a 1d plots of the relative bin errors
	361	//
	362
	363	TH3* hist = GetCorrectionHistogram();
	364	if (!hist)
	365	return 0;
	366
	367	TH1F* target = new TH1F("relerrors", "relerrors", 100, 0, 10);
	368
	369	for (Int_t x=hist->GetXaxis()->FindBin(xmin); x<=hist->GetXaxis()->FindBin(xmax); ++x)
	370	for (Int_t y=hist->GetYaxis()->FindBin(ymin); y<=hist->GetYaxis()->FindBin(ymax); ++y)
	371	for (Int_t z=hist->GetZaxis()->FindBin(zmin); z<=hist->GetZaxis()->FindBin(zmax); ++z)
	372	if (hist->GetBinContent(x, y, z) != 0)
	373	target->Fill(100 * hist->GetBinError(x, y, z) / hist->GetBinContent(x, y, z));
	374
	375	return target;
	376	}
	377