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

/* $Id$ */

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

#include <TH2F.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
  //

  fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);
  fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, Xmin, Xmax, nBinY, Ymin, Ymax);

  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
  //

  fhMeas  = new TH2F(Form("meas_%s",name), Form("meas_%s",title),  nBinX, X, nBinY, Y);
  fhGene  = new TH2F(Form("gene_%s",name), Form("gene_%s",title),  nBinX, X, nBinY, Y);
  fhCorr  = new TH2F(Form("corr_%s",name), Form("corr_%s",title),  nBinX, X, nBinY, Y);

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

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

  // histograms already deleted in base class
}

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

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

//____________________________________________________________________
TH1F* AliCorrectionMatrix2D::Get1DCorrection(Char_t* opt, Float_t min, Float_t max)
{
  //
  // integrate the correction over one variable 
  // 

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

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

    if (min==0 && max==0) {
      meas1D = GetMeasuredHistogram()->ProjectionX();
      gene1D = GetGeneratedHistogram()->ProjectionX();
    }
    else {
      AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including y-bins %d to %d \n", binMin, binMax));

      meas1D = GetMeasuredHistogram()->ProjectionX("pm",binMin,binMax);
      gene1D = GetGeneratedHistogram()->ProjectionX("pg",binMin,binMax);
    }
  }
  if (strcmp(opt,"y")==0) {
    Int_t binMin = GetMeasuredHistogram()->GetXaxis()->FindBin(min);
    Int_t binMax = GetMeasuredHistogram()->GetXaxis()->FindBin(max);

    if (min==0 && max==0) {
      meas1D = GetMeasuredHistogram()->ProjectionY();
      gene1D = GetGeneratedHistogram()->ProjectionY();
    }
    else {
      AliDebug(AliLog::kDebug+1, Form("Getting 1D map. Including x-bins %d to %d \n", binMin, binMax));

      meas1D = GetMeasuredHistogram()->ProjectionY("pm", binMin, binMax);
      gene1D = GetGeneratedHistogram()->ProjectionY("pg", binMin, binMax);
    }
  }
  gene1D->Sumw2();

  gene1D->SetName(Form("corr_1D_%s",fName.Data()));
  gene1D->SetTitle(Form("corr_1D_%s",fName.Data()));

  gene1D->Divide(gene1D, meas1D, 1, 1, "B");
  
  return (TH1F*)gene1D;   
}

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

//____________________________________________________________________
void AliCorrectionMatrix2D::FillGene(Float_t ax, Float_t ay)
{
  // add value to generated histogram
  GetGeneratedHistogram()->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();  

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