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

/* $Id$ */

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

#include <TH2F.h>
#include <TMath.h>

#include <AliLog.h>

#include "AliCorrectionMatrix2D.h"

//____________________________________________________________________
ClassImp(AliCorrectionMatrix2D)

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

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

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

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

  return *this;
}

//____________________________________________________________________
AliCorrectionMatrix2D::AliCorrectionMatrix2D(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) 
  : AliCorrectionMatrix(name, title)
{
  //
  // constructor
  //

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

  fhMeas  = new TH2F("measured",   Form("%s measured", GetTitle()),   nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  fhGene  = new TH2F("generated",  Form("%s generated", GetTitle()),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  fhCorr  = new TH2F("correction", Form("%s correction", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);

  TH1::AddDirectory(oldStatus);

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

//____________________________________________________________________
AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
				       Int_t nBinX, Float_t *X, Int_t nBinY, Float_t *Y) 
  : AliCorrectionMatrix(name, title)
{
  //
  // constructor
  //

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

	fhMeas  = new TH2F("measured",   Form("%s measured",title),   nBinX, X, nBinY, Y);
  fhGene  = new TH2F("generated",  Form("%s generated",title),  nBinX, X, nBinY, Y);
  fhCorr  = new TH2F("correction", Form("%s correction",title), nBinX, X, nBinY, Y);

  TH1::AddDirectory(oldStatus);

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

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

  // histograms already deleted in base class
}

TH2* AliCorrectionMatrix2D::GetGeneratedHistogram() const
{
  // return generated histogram casted to correct type
  return dynamic_cast<TH2*> (fhGene);
}

TH2* AliCorrectionMatrix2D::GetMeasuredHistogram() const
{
  // return measured histogram casted to correct type
  return dynamic_cast<TH2*> (fhMeas);
}

//____________________________________________________________________
TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(const Char_t* opt, Float_t min, Float_t max, Bool_t binomialErrors)
{
  //
  // integrate the correction over one variable 
  // 

  TH1D* meas1D = 0;
  TH1D* gene1D = 0; 

  if (strcmp(opt,"x")==0) {
    Int_t binMin = fhMeas->GetYaxis()->FindBin(min);
    Int_t binMax = fhGene->GetYaxis()->FindBin(max);

    if (min>=max) {
      meas1D = ((TH2F*)fhMeas)->ProjectionX();
      gene1D = ((TH2F*)fhGene)->ProjectionX();
    }
    else {
      Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);

      meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_x_pm", GetName()),binMin,binMax);
      gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_x_pg", GetName()),binMin,binMax);
    }
  }
  else if (strcmp(opt,"y")==0) {
    Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
    Int_t binMax = fhMeas->GetXaxis()->FindBin(max);

    if (min>=max) {
      meas1D = ((TH2F*)fhMeas)->ProjectionY();
      gene1D = ((TH2F*)fhGene)->ProjectionY();
    }
    else {
      Printf("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax);

      meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_y_pm", GetName()), binMin, binMax);
      gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_y_pg", GetName()), binMin, binMax);
    }
  }
  else {
    Printf("ERROR: Invalid option");
    return 0;
  }

  if (!binomialErrors)
  {
    // set the errors on gene manually, and clear the ones on meas.
    gene1D->Sumw2();
    for (Int_t bin=0; bin <= gene1D->GetNbinsX()+1; bin++)
    {
      gene1D->SetBinError(bin, TMath::Sqrt(gene1D->GetBinContent(bin)));
      meas1D->SetBinError(bin, 0);
    }
  }
  
  gene1D->SetName(Form("corr_1D_%s",fName.Data()));
  gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
 
  TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
  divided->Reset();
  
  divided->Divide(gene1D, meas1D, 1, 1, (binomialErrors) ? "B" : "");

  Printf("%p %p", gene1D, meas1D);
  
  return (TH1F*)divided;   
}

//____________________________________________________________________
void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
{
  // add value to measured histogram
  ((TH2F*)fhMeas)->Fill(ax, ay);
}

//____________________________________________________________________
void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
{
  // add value to generated histogram
  ((TH2F*)fhGene)->Fill(ax, ay);
}

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

//____________________________________________________________________
void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
{
  // remove edges of correction histogram by removing 
  // - bins with content less than cut
  // - bins next to bins with zero bin content
  
  Int_t nBinsX = fhCorr->GetNbinsX();
  Int_t nBinsY = fhCorr->GetNbinsY();

  // set bin content to zero for bins with content smaller cut
  for (Int_t bx=0; bx<=nBinsX; bx++) {
    for (Int_t by=0; by<=nBinsY; by++) {
      if (fhCorr->GetBinContent(bx,by)>cut) {
	  fhCorr->SetBinContent(bx,by,0);
	  fhCorr->SetBinError(bx,by,0);
      }
    }
  }

  // set bin content to zero for bins next to bins with zero
  TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
  tmp->Reset();
  
  Bool_t done = kFALSE;
  Int_t nBinsXCount = 0;
  Int_t nBinsYCount = 0;
  while (!done) {    
    if (nBinsXCount<nBinsXedge) 
      for (Int_t bx=0; bx<=nBinsX; bx++) {
	for (Int_t by=0; by<=nBinsY; by++) {
	  if ((fhCorr->GetBinContent(bx+1,by)==0)|| 
	      (fhCorr->GetBinContent(bx-1,by)==0))
	    tmp->SetBinContent(bx,by,1);	
	  
	}
      }
    if (nBinsYCount<nBinsYedge) 
      for (Int_t bx=0; bx<=nBinsX; bx++) {
	for (Int_t by=0; by<=nBinsY; by++) {
	  if ((fhCorr->GetBinContent(bx,by+1)==0)|| 
	      (fhCorr->GetBinContent(bx,by-1)==0))
	    tmp->SetBinContent(bx,by,1);	
	}
      }    
    for (Int_t bx=0; bx<=nBinsX; bx++) {
      for (Int_t by=0; by<=nBinsY; by++) {
	if (tmp->GetBinContent(bx,by)==1) {
	  fhCorr->SetBinContent(bx,by,0);
	  fhCorr->SetBinError(bx,by,0);
	}
      }
    }
    nBinsXCount++;
    nBinsYCount++;
    if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
  }
  tmp->Delete();  

}

//____________________________________________________________________
void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
{
	// rebins the histograms, recalculates the correction

        GetGeneratedHistogram()->Rebin2D(x, y);
        GetMeasuredHistogram()->Rebin2D(x, y);
        GetCorrectionHistogram()->Rebin2D(x, y);
        Divide();
}
Commit	Line	Data
bf21645b	1	/* $Id$ */
	2
	3	// ------------------------------------------------------
	4	//
	5	// Class to handle 2d-corrections.
	6	//
	7	// ------------------------------------------------------
	8	//
	9
	10	#include <TH2F.h>
69b09e3b	11	#include <TMath.h>
bf21645b	12
	13	#include <AliLog.h>
	14
	15	#include "AliCorrectionMatrix2D.h"
	16
	17	//____________________________________________________________________
	18	ClassImp(AliCorrectionMatrix2D)
	19
	20	//____________________________________________________________________
	21	AliCorrectionMatrix2D::AliCorrectionMatrix2D() :
	22	AliCorrectionMatrix()
	23	{
	24	// default constructor
	25	}
	26
	27	//____________________________________________________________________
	28	AliCorrectionMatrix2D::AliCorrectionMatrix2D(const AliCorrectionMatrix2D& c)
	29	: AliCorrectionMatrix(c)
	30	{
	31	// copy constructor
	32	((AliCorrectionMatrix2D &)c).Copy(*this);
	33	}
	34
61385583	35	//____________________________________________________________________
	36	AliCorrectionMatrix2D &AliCorrectionMatrix2D::operator=(const AliCorrectionMatrix2D &c)
	37	{
	38	// assigment operator
	39
	40	if (this != &c)
	41	((AliCorrectionMatrix2D &) c).Copy(*this);
	42
	43	return *this;
	44	}
	45
bf21645b	46	//____________________________________________________________________
	47	AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
	48	Int_t nBinX, Float_t Xmin, Float_t Xmax,
	49	Int_t nBinY, Float_t Ymin, Float_t Ymax)
	50	: AliCorrectionMatrix(name, title)
	51	{
	52	//
	53	// constructor
	54	//
	55
29771dc8	56	// do not add this hists to the directory
	57	Bool_t oldStatus = TH1::AddDirectoryStatus();
	58	TH1::AddDirectory(kFALSE);
	59
82ac9fc4	60	fhMeas = new TH2F("measured", Form("%s measured", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
	61	fhGene = new TH2F("generated", Form("%s generated", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
	62	fhCorr = new TH2F("correction", Form("%s correction", GetTitle()), nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
29771dc8	63
29771dc8	64	TH1::AddDirectory(oldStatus);
bf21645b	65
	66	fhMeas->Sumw2();
	67	fhGene->Sumw2();
	68	fhCorr->Sumw2();
	69	}
	70
	71	//____________________________________________________________________
	72	AliCorrectionMatrix2D::AliCorrectionMatrix2D(const Char_t* name, const Char_t* title,
	73	Int_t nBinX, Float_t X, Int_t nBinY, Float_t Y)
	74	: AliCorrectionMatrix(name, title)
	75	{
	76	//
	77	// constructor
	78	//
	79
29771dc8	80	// do not add this hists to the directory
	81	Bool_t oldStatus = TH1::AddDirectoryStatus();
	82	TH1::AddDirectory(kFALSE);
bf21645b	83
29771dc8	84	fhMeas = new TH2F("measured", Form("%s measured",title), nBinX, X, nBinY, Y);
	85	fhGene = new TH2F("generated", Form("%s generated",title), nBinX, X, nBinY, Y);
	86	fhCorr = new TH2F("correction", Form("%s correction",title), nBinX, X, nBinY, Y);
	87
	88	TH1::AddDirectory(oldStatus);
	89
	90	fhMeas->Sumw2();
bf21645b	91	fhGene->Sumw2();
	92	fhCorr->Sumw2();
	93	}
	94
	95	//____________________________________________________________________
	96	AliCorrectionMatrix2D::~AliCorrectionMatrix2D()
	97	{
	98	//
	99	// destructor
	100	//
	101
	102	// histograms already deleted in base class
	103	}
	104
5a6310fe	105	TH2* AliCorrectionMatrix2D::GetGeneratedHistogram() const
bf21645b	106	{
bf21645b	107	// return generated histogram casted to correct type
5a6310fe	108	return dynamic_cast<TH2*> (fhGene);
bf21645b	109	}
bf21645b	110
5a6310fe	111	TH2* AliCorrectionMatrix2D::GetMeasuredHistogram() const
bf21645b	112	{
bf21645b	113	// return measured histogram casted to correct type
5a6310fe	114	return dynamic_cast<TH2*> (fhMeas);
bf21645b	115	}
	116
	117	//____________________________________________________________________
69b09e3b	118	TH1* AliCorrectionMatrix2D::Get1DCorrectionHistogram(const Char_t* opt, Float_t min, Float_t max, Bool_t binomialErrors)
bf21645b	119	{
	120	//
	121	// integrate the correction over one variable
	122	//
	123
	124	TH1D* meas1D = 0;
	125	TH1D* gene1D = 0;
	126
	127	if (strcmp(opt,"x")==0) {
06e4b91b	128	Int_t binMin = fhMeas->GetYaxis()->FindBin(min);
06e4b91b	129	Int_t binMax = fhGene->GetYaxis()->FindBin(max);
7584d357	130
06e4b91b	131	if (min>=max) {
	132	meas1D = ((TH2F*)fhMeas)->ProjectionX();
	133	gene1D = ((TH2F*)fhGene)->ProjectionX();
7584d357	134	}
7584d357	135	else {
dd367a14	136	Printf("Getting 1D map. Including y-bins %d to %d", binMin, binMax);
7584d357	137
5a6310fe	138	meas1D = ((TH2F*)fhMeas)->ProjectionX(Form("%s_x_pm", GetName()),binMin,binMax);
5a6310fe	139	gene1D = ((TH2F*)fhGene)->ProjectionX(Form("%s_x_pg", GetName()),binMin,binMax);
7584d357	140	}
bf21645b	141	}
5a6310fe	142	else if (strcmp(opt,"y")==0) {
06e4b91b	143	Int_t binMin = fhMeas->GetXaxis()->FindBin(min);
06e4b91b	144	Int_t binMax = fhMeas->GetXaxis()->FindBin(max);
dd367a14	145
06e4b91b	146	if (min>=max) {
	147	meas1D = ((TH2F*)fhMeas)->ProjectionY();
	148	gene1D = ((TH2F*)fhGene)->ProjectionY();
7584d357	149	}
7584d357	150	else {
5a6310fe	151	Printf("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax);
dd367a14	152
5a6310fe	153	meas1D = ((TH2F*)fhMeas)->ProjectionY(Form("%s_y_pm", GetName()), binMin, binMax);
5a6310fe	154	gene1D = ((TH2F*)fhGene)->ProjectionY(Form("%s_y_pg", GetName()), binMin, binMax);
7584d357	155	}
bf21645b	156	}
5a6310fe	157	else {
	158	Printf("ERROR: Invalid option");
	159	return 0;
	160	}
	161
69b09e3b	162	if (!binomialErrors)
	163	{
	164	// set the errors on gene manually, and clear the ones on meas.
	165	gene1D->Sumw2();
	166	for (Int_t bin=0; bin <= gene1D->GetNbinsX()+1; bin++)
	167	{
	168	gene1D->SetBinError(bin, TMath::Sqrt(gene1D->GetBinContent(bin)));
	169	meas1D->SetBinError(bin, 0);
	170	}
	171	}
	172
bf21645b	173	gene1D->SetName(Form("corr_1D_%s",fName.Data()));
bf21645b	174	gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));
51f6de65	175
	176	TH1* divided = (TH1*) gene1D->Clone(Form("corr_1D_%s",fName.Data()));
	177	divided->Reset();
69b09e3b	178
69b09e3b	179	divided->Divide(gene1D, meas1D, 1, 1, (binomialErrors) ? "B" : "");
bf21645b	180
51f6de65	181	Printf("%p %p", gene1D, meas1D);
bf21645b	182
51f6de65	183	return (TH1F*)divided;
bf21645b	184	}
	185
	186	//____________________________________________________________________
	187	void AliCorrectionMatrix2D::FillMeas(Float_t ax, Float_t ay)
	188	{
	189	// add value to measured histogram
06e4b91b	190	((TH2F*)fhMeas)->Fill(ax, ay);
bf21645b	191	}
	192
	193	//____________________________________________________________________
	194	void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
	195	{
	196	// add value to generated histogram
06e4b91b	197	((TH2F*)fhGene)->Fill(ax, ay);
bf21645b	198	}
	199
	200	//____________________________________________________________________
	201	Float_t AliCorrectionMatrix2D::GetCorrection(Float_t ax, Float_t ay) const
	202	{
	203	// returns a value of the correction map
	204	return fhCorr->GetBinContent(fhCorr->FindBin(ax,ay));
	205	}
	206
	207	//____________________________________________________________________
	208	void AliCorrectionMatrix2D::RemoveEdges(Float_t cut, Int_t nBinsXedge, Int_t nBinsYedge)
	209	{
	210	// remove edges of correction histogram by removing
	211	// - bins with content less than cut
	212	// - bins next to bins with zero bin content
	213
	214	Int_t nBinsX = fhCorr->GetNbinsX();
	215	Int_t nBinsY = fhCorr->GetNbinsY();
	216
	217	// set bin content to zero for bins with content smaller cut
	218	for (Int_t bx=0; bx<=nBinsX; bx++) {
	219	for (Int_t by=0; by<=nBinsY; by++) {
	220	if (fhCorr->GetBinContent(bx,by)>cut) {
	221	fhCorr->SetBinContent(bx,by,0);
	222	fhCorr->SetBinError(bx,by,0);
	223	}
	224	}
	225	}
	226
	227	// set bin content to zero for bins next to bins with zero
	228	TH2F* tmp = (TH2F*)fhCorr->Clone("tmp");
	229	tmp->Reset();
	230
	231	Bool_t done = kFALSE;
	232	Int_t nBinsXCount = 0;
	233	Int_t nBinsYCount = 0;
	234	while (!done) {
	235	if (nBinsXCount<nBinsXedge)
	236	for (Int_t bx=0; bx<=nBinsX; bx++) {
	237	for (Int_t by=0; by<=nBinsY; by++) {
	238	if ((fhCorr->GetBinContent(bx+1,by)==0)\|\|
	239	(fhCorr->GetBinContent(bx-1,by)==0))
	240	tmp->SetBinContent(bx,by,1);
	241
	242	}
	243	}
	244	if (nBinsYCount<nBinsYedge)
	245	for (Int_t bx=0; bx<=nBinsX; bx++) {
	246	for (Int_t by=0; by<=nBinsY; by++) {
	247	if ((fhCorr->GetBinContent(bx,by+1)==0)\|\|
	248	(fhCorr->GetBinContent(bx,by-1)==0))
	249	tmp->SetBinContent(bx,by,1);
	250	}
	251	}
	252	for (Int_t bx=0; bx<=nBinsX; bx++) {
	253	for (Int_t by=0; by<=nBinsY; by++) {
	254	if (tmp->GetBinContent(bx,by)==1) {
	255	fhCorr->SetBinContent(bx,by,0);
	256	fhCorr->SetBinError(bx,by,0);
	257	}
	258	}
	259	}
	260	nBinsXCount++;
	261	nBinsYCount++;
262	if ((nBinsXCount>=nBinsXedge)&&(nBinsYCount>=nBinsYedge)) done=kTRUE;
263	}
264	tmp->Delete();
265
266	}
267
745d6088	268	//____________________________________________________________________
	269	void AliCorrectionMatrix2D::Rebin(Int_t x, Int_t y)
	270	{
	271	// rebins the histograms, recalculates the correction
	272
	273	GetGeneratedHistogram()->Rebin2D(x, y);
	274	GetMeasuredHistogram()->Rebin2D(x, y);
	275	GetCorrectionHistogram()->Rebin2D(x, y);
	276	Divide();
	277	}