[u/mrichter/AliRoot.git] / TPC / AliH2F.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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.1.4.2  2000/04/10 11:32:37  kowal2

"ROOT"-based class with some extra functionality

*/

//----------------------------------------------------------------------------
//  Author:   Marian Ivanov
//
//  Implementation of class AliH2F
//
//-----------------------------------------------------------------------------

#include "AliH2F.h"
#include "TClonesArray.h"
#include "AliTPC.h"
#include "TRandom.h"
#include "AliCluster.h"
#include "AliClusterFinder.h"
//*KEEP,TMath.

// other include files follow here


ClassImp(AliH2F)
//***********************************************************************
//***********************************************************************
//***********************************************************************
//***********************************************************************
AliH2F::AliH2F():TH2F() 
{
  fSigmaX2 = 1.;
  fSigmaY2 = 1.;
  fdx = 1;
  fdy=1;
  fFitMatrix.ResizeTo(5,1);  
}
AliH2F::AliH2F(const Text_t *name,const Text_t *title,
		       Int_t nbinsx,Axis_t xlow,Axis_t xup
		       ,Int_t nbinsy,Axis_t ylow,Axis_t yup):
  TH2F(name,title,nbinsx,xlow,xup
       ,nbinsy,ylow,yup)
{
  fSigmaX2 = 1.;
  fSigmaY2 = 1.;
  fdx = 1;
  fdy=1; 
  fFitMatrix.ResizeTo(5,1);  
}
     
AliH2F::~AliH2F() 
{
}

AliH2F::AliH2F(const AliH2F &) 
{
}

AliH2F & AliH2F::operator = (const AliH2F &) 
{
   return *this;
}

TClonesArray * AliH2F::FindPeaks(Float_t threshold, Float_t noise)
{
  //find peaks and write it in form of AliTPCcluster to array
    
  //firstly we need to create object for cluster finding
  //and fill it with contents of histogram
  AliClusterFinder cfinder;
  cfinder.SetThreshold(threshold);
  cfinder.SetNoise(noise);
  cfinder.GetHisto(this);
  return cfinder.FindPeaks3();
}

void AliH2F::ClearSpectrum()
{
  Int_t dimx =  fXaxis.GetNbins();
  Int_t dimy =  fYaxis.GetNbins();
  for (Int_t i = 0 ;i<dimx;i++)
    for (Int_t j = 0 ;j<dimy;j++) 
      {
	SetCellContent(i,j,0);
	SetCellError(i,j,0);
      }
}


void AliH2F::AddNoise(Float_t sn)
{
  Int_t dimx =  fXaxis.GetNbins();
  Int_t dimy =  fYaxis.GetNbins();
  for (Int_t i = 0 ;i<dimx;i++)
    for (Int_t j = 0 ;j<dimy;j++) 
      {
        Float_t noise = gRandom->Gaus(0,sn);
	Float_t oldv  =GetCellContent(i,j);
	Float_t olds  =GetCellError(i,j);
	if (noise >0)
	  {
	    SetCellContent(i,j,noise+oldv);
	    SetCellError(i,j,TMath::Sqrt((noise*noise+olds*olds)));
	  }
      }
}
void AliH2F::AddGauss(Float_t x, Float_t y, 
			  Float_t sx, Float_t sy, Float_t max)
{
  //transform to histogram coordinata  
  Int_t dimx =  fXaxis.GetNbins();
  Int_t dimy =  fYaxis.GetNbins();
  Float_t dx =(GetXaxis()->GetXmax()-GetXaxis()->GetXmin())/Float_t(dimx);
  Float_t dy =(GetYaxis()->GetXmax()-GetYaxis()->GetXmin())/Float_t(dimy);  
  //  x=(x-GetXaxis()->GetXmin())/dx;
  //y=(y-GetYaxis()->GetXmin())/dy;
  sx/=dx;
  sy/=dy;

  
  for (Int_t i = 0 ;i<dimx;i++)
    for (Int_t j = 0 ;j<dimy;j++) 
      {
	Float_t x2 =GetXaxis()->GetBinCenter(i+1);
	Float_t y2 =GetYaxis()->GetBinCenter(j+1);
	Float_t dx2 = (x2-x)*(x2-x);
        Float_t dy2 = (y2-y)*(y2-y);
        Float_t amp =max*exp(-(dx2/(2*sx*sx)+dy2/(2*sy*sy)));
	//Float_t oldv  =GetCellContent(i+1,j+1);
	//	SetCellContent(i+1,j+1,amp+oldv);
	Fill(x2,y2,amp);
      }
}

void AliH2F::ClearUnderTh(Int_t threshold)
{
  Int_t dimx =  fXaxis.GetNbins();
  Int_t dimy =  fYaxis.GetNbins();
  for (Int_t i = 0 ;i<=dimx;i++)
    for (Int_t j = 0 ;j<=dimy;j++) 
      {	
	Float_t oldv  =GetCellContent(i,j);
        if (oldv <threshold)
	  SetCellContent(i,j,0);
      }
}

void AliH2F::Round()
{
  Int_t dimx =  fXaxis.GetNbins();
  Int_t dimy =  fYaxis.GetNbins();
  for (Int_t i = 0 ;i<=dimx;i++)
    for (Int_t j = 0 ;j<=dimy;j++) 
      {	
	Float_t oldv  =GetCellContent(i,j);
        oldv=(Int_t)oldv;
	SetCellContent(i,j,oldv);
      }
}


AliH2F *AliH2F::GetSubrange2d(Float_t xmin, Float_t xmax, 
				      Float_t ymin, Float_t ymax)
{
  //this function return pointer to the new created 
  //histogram which is subhistogram of the 
  //calculate number
  //subhistogram range must be inside histogram

  if (xmax<=xmin) {
    xmin=fXaxis.GetXmin();
    xmax=fXaxis.GetXmax();
  }
  if (ymax<=ymin) {
     ymin=fYaxis.GetXmin();
     ymax=fYaxis.GetXmax();
  }

  Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  * 
		   Float_t(fXaxis.GetNbins()));
  Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  * 
		   Float_t(fYaxis.GetNbins()));
  TString  t1 = fName ;
  TString  t2 = fTitle ;
  t1+="_subrange";
  t2+="_subrange";
  const  Text_t * tt1 = t1;
  const Text_t * tt2 = t2;
  
  AliH2F * sub = new AliH2F(tt1,tt2,nx,xmin,xmax,ny,ymin,ymax); 
  
  Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
		    (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
  Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
		    (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
  for (Int_t i=0;i<nx;i++)
    for (Int_t j=0;j<ny;j++)
      {
	Int_t index1 = GetBin(i1+i,i2+j);
	//        Int_t index2 = sub->GetBin(i,j);
        Float_t val = GetBinContent(index1);
	//        sub->SetBinContent(index2,val);
	//        Float_t err = GetBinError(index1);
        //sub->SetBinError(index2,GetBinError(index1));
        sub->SetCellContent(i,j,val);
      }  
   return sub;
}

TH1F *AliH2F::GetAmplitudes(Float_t zmin, Float_t zmax, Float_t th, Float_t xmin, Float_t xmax, 
				      Float_t ymin, Float_t ymax)
{
  //this function return pointer to the new created 
  //histogram which is subhistogram of the 
  //calculate number
  //subhistogram range must be inside histogram
 
  if (xmax<=xmin) {
    xmin=fXaxis.GetXmin();
    xmax=fXaxis.GetXmax();
  }
  if (ymax<=ymin) {
     ymin=fYaxis.GetXmin();
     ymax=fYaxis.GetXmax();
  }
  Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  * 
		   Float_t(fXaxis.GetNbins()));
  Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  * 
		   Float_t(fYaxis.GetNbins()));
  TString  t1 = fName ;
  TString  t2 = fTitle ;
  t1+="_amplitudes";
  t2+="_amplitudes";
  const  Text_t * tt1 = t1;
  const Text_t * tt2 = t2;
  
  TH1F * h = new TH1F(tt1,tt2,100,zmin,zmax); 
  
  Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
		    (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
  Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
		    (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
  for (Int_t i=0;i<nx;i++)
    for (Int_t j=0;j<ny;j++)
      {
	Int_t index1 = GetBin(i1+i,i2+j);
        Float_t val = GetBinContent(index1);
        if (val>th) h->Fill(val);
      }  
   return h;
}

Float_t   AliH2F::GetOccupancy(Float_t th , Float_t xmin, Float_t xmax, 
			     Float_t ymin, Float_t ymax)
{
  //this function return pointer to the new created 
  //histogram which is subhistogram of the 
  //calculate number
  //subhistogram range must be inside histogram
 
  if (xmax<=xmin) {
    xmin=fXaxis.GetXmin();
    xmax=fXaxis.GetXmax();
  }
  if (ymax<=ymin) {
     ymin=fYaxis.GetXmin();
     ymax=fYaxis.GetXmax();
  }
  Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  * 
		   Float_t(fXaxis.GetNbins()));
  Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  * 
		   Float_t(fYaxis.GetNbins()));
 
  Int_t over =0; 
  Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
		    (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
  Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
		    (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
  for (Int_t i=0;i<nx;i++)
    for (Int_t j=0;j<ny;j++)
      {
	Int_t index1 = GetBin(i1+i,i2+j);
        Float_t val = GetBinContent(index1);
        if (val>th) over++;
      }  
  Int_t  all = nx*ny;
  if (all>0)  return Float_t(over)/Float_t(all);
  else 
    return 0;
}

//TH1F * AliH2F::GetSubrange1dx(Float_t xmin, Float_t xmax, Float_t y)
//{
//
//}
//TH1F * AliH2F::GetSubrange1dy(Float_t x, Float_t ymin, Float_t ymax)
//{
//
//}

void AliH2F::SetSmoothSigma(Float_t sigmaX, Float_t sigmaY)
{
  Float_t wx = (fXaxis.GetXmax()-fXaxis.GetXmin()) / Float_t(fXaxis.GetNbins());
  Float_t wx2 = wx*wx;
  Float_t wy = (fYaxis.GetXmax()-fYaxis.GetXmin()) / Float_t(fYaxis.GetNbins()) ;
  Float_t wy2 =wy*wy;
  fSigmaX2 = sigmaX*sigmaX/wx2;
  fSigmaY2 = sigmaY*sigmaX/wy2;  
}

void AliH2F::SetSmoothRange(Float_t dx,Float_t dy)
{
  Float_t wx = (fXaxis.GetXmax()-fXaxis.GetXmin()) / Float_t(fXaxis.GetNbins());
  Float_t wy = (fYaxis.GetXmax()-fYaxis.GetXmin()) / Float_t(fYaxis.GetNbins());
  fdx = Int_t(dx/wx);
  fdy = Int_t(dy/wy);
}

TMatrix *  AliH2F::GetFitMatrix(Int_t irow, Int_t icolumn) 
{
  SmoothCell(irow,icolumn);
  return &fFitMatrix;
}

void AliH2F::SmoothCell(Int_t irow, Int_t icolumn)
{
  //calculate interpolation around point irow icollumn
  //fitting by surface of second oreder
  Int_t index;  
  Float_t x,y,z,x2,x3,x4,y2,y3,y4;
  Float_t sz = 0.0,szx = 0.0 ,szy = 0.0 ,szx2 = 0.0,szy2 = 0.0;
  Float_t sx = 0.,sx2 = 0.,sx3 = 0.,sx4 = 0.;
  Float_t sy = 0.,sy2 = 0. ,sy3 = 0. ,sy4 = 0.;
  Float_t sxy = 0.,sxy2 = 0. ,sx2y =0.0 ,sx2y2 =0.0;
  Float_t w;
  Float_t sumweight = 0;
  for (Int_t i = -fdx; i <=fdx; i++)  
    for (Int_t j = -fdy; j <=fdy; j++)  
      {
	index = GetBin(irow+i,icolumn+j);
        w = 1/(Float_t(i*i)+fSigmaX2)*1/(Float_t(j*j)+fSigmaY2);
        z = GetBinContent(index);
        x = irow+i;
        x2 = x*x;           x3 = x2*x;         x4 = x2*x2;
        y = icolumn+j;
        y2   = y*y;         y3   = y2*y;       y4   = y2*y2; 
	sz  += z*w;         sx  += x*w;        sy  += y*w;
        szx += z*x*w;       szy += z*y*w;      szx2+= z*x2*w;     szy2 += z*y3*w;
        sx2 += x2*w;        sx3 += x3*w;       sx4 += x4*w;
        sy2 += y2*w;        sy3 += y3*w;       sy4 += y4*w;
        sxy += x*y*w;       sxy2+= x*y2*w;     sx2y+=x2*y*w;      sx2y2+= x2*y2*w;
        sumweight +=w;
      }
  TMatrix mat(5,5);
  if (!fFitMatrix.IsValid()) fFitMatrix.ResizeTo(5,1);
  
  fFitMatrix(0,0) = sz;  
  fFitMatrix(1,0) = szx;
  fFitMatrix(2,0) = szy;
  fFitMatrix(3,0) = szx2;
  fFitMatrix(4,0) = szy2;
  
  mat(0,0) = sumweight;
                  mat(0,1) = sx;    mat(0,2) = sy;  mat(0,3) = sx2;  mat(0,4) =sy2; 
  mat(1,0) = sx;  mat(1,1) = sx2;   mat(1,2) = sxy; mat(1,3) = sx3;  mat(1,4) =sxy2;
  mat(2,0) = sy;  mat(2,1) = sxy;   mat(2,2) = sy2; mat(2,3) = sx2y; mat(2,4) =sy3 ;
  mat(3,0) = sx2; mat(3,1) = sx3;   mat(3,2) = sx2y;mat(3,3) = sx4 ; mat(3,4) =sx2y2;
  mat(4,0) = sy2; mat(4,1) = sxy2;  mat(4,2) = sy3; mat(4,3) = sx2y2;mat(4,4) =sy4;
}
Commit	Line	Data
cc80f89e	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	/*
	17	$Log$
	18	Revision 1.1.4.2 2000/04/10 11:32:37 kowal2
	19
	20	"ROOT"-based class with some extra functionality
	21
	22	*/
	23
	24	//----------------------------------------------------------------------------
	25	// Author: Marian Ivanov
	26	//
	27	// Implementation of class AliH2F
	28	//
	29	//-----------------------------------------------------------------------------
	30
	31	#include "AliH2F.h"
	32	#include "TClonesArray.h"
	33	#include "AliTPC.h"
	34	#include "TRandom.h"
	35	#include "AliCluster.h"
	36	#include "AliClusterFinder.h"
	37	//*KEEP,TMath.
	38
	39	// other include files follow here
	40
	41
	42	ClassImp(AliH2F)
	43	//***********************************************************************
	44	//***********************************************************************
	45	//***********************************************************************
	46	//***********************************************************************
	47	AliH2F::AliH2F():TH2F()
	48	{
	49	fSigmaX2 = 1.;
	50	fSigmaY2 = 1.;
	51	fdx = 1;
	52	fdy=1;
	53	fFitMatrix.ResizeTo(5,1);
	54	}
	55	AliH2F::AliH2F(const Text_t name,const Text_t title,
	56	Int_t nbinsx,Axis_t xlow,Axis_t xup
	57	,Int_t nbinsy,Axis_t ylow,Axis_t yup):
	58	TH2F(name,title,nbinsx,xlow,xup
	59	,nbinsy,ylow,yup)
	60	{
	61	fSigmaX2 = 1.;
	62	fSigmaY2 = 1.;
	63	fdx = 1;
	64	fdy=1;
65	fFitMatrix.ResizeTo(5,1);
66	}
67
68	AliH2F::~AliH2F()
69	{
70	}
71
72	AliH2F::AliH2F(const AliH2F &)
73	{
74	}
75
76	AliH2F & AliH2F::operator = (const AliH2F &)
77	{
78	return *this;
79	}
80
81	TClonesArray * AliH2F::FindPeaks(Float_t threshold, Float_t noise)
82	{
83	//find peaks and write it in form of AliTPCcluster to array
84
85	//firstly we need to create object for cluster finding
86	//and fill it with contents of histogram
87	AliClusterFinder cfinder;
88	cfinder.SetThreshold(threshold);
89	cfinder.SetNoise(noise);
90	cfinder.GetHisto(this);
91	return cfinder.FindPeaks3();
92	}
93
94	void AliH2F::ClearSpectrum()
95	{
96	Int_t dimx = fXaxis.GetNbins();
97	Int_t dimy = fYaxis.GetNbins();
98	for (Int_t i = 0 ;i<dimx;i++)
99	for (Int_t j = 0 ;j<dimy;j++)
100	{
101	SetCellContent(i,j,0);
102	SetCellError(i,j,0);
103	}
104	}
105
106
107	void AliH2F::AddNoise(Float_t sn)
108	{
109	Int_t dimx = fXaxis.GetNbins();
110	Int_t dimy = fYaxis.GetNbins();
111	for (Int_t i = 0 ;i<dimx;i++)
112	for (Int_t j = 0 ;j<dimy;j++)
113	{
114	Float_t noise = gRandom->Gaus(0,sn);
115	Float_t oldv =GetCellContent(i,j);
116	Float_t olds =GetCellError(i,j);
117	if (noise >0)
118	{
119	SetCellContent(i,j,noise+oldv);
120	SetCellError(i,j,TMath::Sqrt((noisenoise+oldsolds)));
121	}
122	}
123	}
124	void AliH2F::AddGauss(Float_t x, Float_t y,
125	Float_t sx, Float_t sy, Float_t max)
126	{
127	//transform to histogram coordinata
128	Int_t dimx = fXaxis.GetNbins();
129	Int_t dimy = fYaxis.GetNbins();
130	Float_t dx =(GetXaxis()->GetXmax()-GetXaxis()->GetXmin())/Float_t(dimx);
131	Float_t dy =(GetYaxis()->GetXmax()-GetYaxis()->GetXmin())/Float_t(dimy);
132	// x=(x-GetXaxis()->GetXmin())/dx;
133	//y=(y-GetYaxis()->GetXmin())/dy;
134	sx/=dx;
135	sy/=dy;
136
137
138	for (Int_t i = 0 ;i<dimx;i++)
139	for (Int_t j = 0 ;j<dimy;j++)
140	{
141	Float_t x2 =GetXaxis()->GetBinCenter(i+1);
142	Float_t y2 =GetYaxis()->GetBinCenter(j+1);
143	Float_t dx2 = (x2-x)*(x2-x);
144	Float_t dy2 = (y2-y)*(y2-y);
145	Float_t amp =maxexp(-(dx2/(2sxsx)+dy2/(2sy*sy)));
146	//Float_t oldv =GetCellContent(i+1,j+1);
147	// SetCellContent(i+1,j+1,amp+oldv);
148	Fill(x2,y2,amp);
149	}
150	}
151
152	void AliH2F::ClearUnderTh(Int_t threshold)
153	{
154	Int_t dimx = fXaxis.GetNbins();
155	Int_t dimy = fYaxis.GetNbins();
156	for (Int_t i = 0 ;i<=dimx;i++)
157	for (Int_t j = 0 ;j<=dimy;j++)
158	{
159	Float_t oldv =GetCellContent(i,j);
160	if (oldv <threshold)
161	SetCellContent(i,j,0);
162	}
163	}
164
165	void AliH2F::Round()
166	{
167	Int_t dimx = fXaxis.GetNbins();
168	Int_t dimy = fYaxis.GetNbins();
169	for (Int_t i = 0 ;i<=dimx;i++)
170	for (Int_t j = 0 ;j<=dimy;j++)
171	{
172	Float_t oldv =GetCellContent(i,j);
173	oldv=(Int_t)oldv;
174	SetCellContent(i,j,oldv);
175	}
176	}
177
178
179
180	AliH2F *AliH2F::GetSubrange2d(Float_t xmin, Float_t xmax,
181	Float_t ymin, Float_t ymax)
182	{
183	//this function return pointer to the new created
184	//histogram which is subhistogram of the
185	//calculate number
186	//subhistogram range must be inside histogram
187
188	if (xmax<=xmin) {
189	xmin=fXaxis.GetXmin();
190	xmax=fXaxis.GetXmax();
191	}
192	if (ymax<=ymin) {
193	ymin=fYaxis.GetXmin();
194	ymax=fYaxis.GetXmax();
195	}
196
197	Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin()) *
198	Float_t(fXaxis.GetNbins()));
199	Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin()) *
200	Float_t(fYaxis.GetNbins()));
201	TString t1 = fName ;
202	TString t2 = fTitle ;
203	t1+="_subrange";
204	t2+="_subrange";
205	const Text_t * tt1 = t1;
206	const Text_t * tt2 = t2;
207
208	AliH2F * sub = new AliH2F(tt1,tt2,nx,xmin,xmax,ny,ymin,ymax);
209
210	Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
211	(fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
212	Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
213	(fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
214	for (Int_t i=0;i<nx;i++)
215	for (Int_t j=0;j<ny;j++)
216	{
217	Int_t index1 = GetBin(i1+i,i2+j);
218	// Int_t index2 = sub->GetBin(i,j);
219	Float_t val = GetBinContent(index1);
220	// sub->SetBinContent(index2,val);
221	// Float_t err = GetBinError(index1);
222	//sub->SetBinError(index2,GetBinError(index1));
223	sub->SetCellContent(i,j,val);
224	}
225	return sub;
226	}
227
228	TH1F *AliH2F::GetAmplitudes(Float_t zmin, Float_t zmax, Float_t th, Float_t xmin, Float_t xmax,
229	Float_t ymin, Float_t ymax)
230	{
231	//this function return pointer to the new created
232	//histogram which is subhistogram of the
233	//calculate number
234	//subhistogram range must be inside histogram
235
236	if (xmax<=xmin) {
237	xmin=fXaxis.GetXmin();
238	xmax=fXaxis.GetXmax();
239	}
240	if (ymax<=ymin) {
241	ymin=fYaxis.GetXmin();
242	ymax=fYaxis.GetXmax();
243	}
244	Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin()) *
245	Float_t(fXaxis.GetNbins()));
246	Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin()) *
247	Float_t(fYaxis.GetNbins()));
248	TString t1 = fName ;
249	TString t2 = fTitle ;
250	t1+="_amplitudes";
251	t2+="_amplitudes";
252	const Text_t * tt1 = t1;
253	const Text_t * tt2 = t2;
254
255	TH1F * h = new TH1F(tt1,tt2,100,zmin,zmax);
256
257	Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
258	(fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
259	Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
260	(fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
261	for (Int_t i=0;i<nx;i++)
262	for (Int_t j=0;j<ny;j++)
263	{
264	Int_t index1 = GetBin(i1+i,i2+j);
265	Float_t val = GetBinContent(index1);
266	if (val>th) h->Fill(val);
267	}
268	return h;
269	}
270
271	Float_t AliH2F::GetOccupancy(Float_t th , Float_t xmin, Float_t xmax,
272	Float_t ymin, Float_t ymax)
273	{
274	//this function return pointer to the new created
275	//histogram which is subhistogram of the
276	//calculate number
277	//subhistogram range must be inside histogram
278
279	if (xmax<=xmin) {
280	xmin=fXaxis.GetXmin();
281	xmax=fXaxis.GetXmax();
282	}
283	if (ymax<=ymin) {
284	ymin=fYaxis.GetXmin();
285	ymax=fYaxis.GetXmax();
286	}
287	Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin()) *
288	Float_t(fXaxis.GetNbins()));
289	Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin()) *
290	Float_t(fYaxis.GetNbins()));
291
292	Int_t over =0;
293	Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
294	(fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
295	Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
296	(fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
297	for (Int_t i=0;i<nx;i++)
298	for (Int_t j=0;j<ny;j++)
299	{
300	Int_t index1 = GetBin(i1+i,i2+j);
301	Float_t val = GetBinContent(index1);
302	if (val>th) over++;
303	}
304	Int_t all = nx*ny;
305	if (all>0) return Float_t(over)/Float_t(all);
306	else
307	return 0;
308	}
309
310	//TH1F * AliH2F::GetSubrange1dx(Float_t xmin, Float_t xmax, Float_t y)
311	//{
312	//
313	//}
314	//TH1F * AliH2F::GetSubrange1dy(Float_t x, Float_t ymin, Float_t ymax)
315	//{
316	//
317	//}
318
319	void AliH2F::SetSmoothSigma(Float_t sigmaX, Float_t sigmaY)
320	{
321	Float_t wx = (fXaxis.GetXmax()-fXaxis.GetXmin()) / Float_t(fXaxis.GetNbins());
322	Float_t wx2 = wx*wx;
323	Float_t wy = (fYaxis.GetXmax()-fYaxis.GetXmin()) / Float_t(fYaxis.GetNbins()) ;
324	Float_t wy2 =wy*wy;
325	fSigmaX2 = sigmaX*sigmaX/wx2;
326	fSigmaY2 = sigmaY*sigmaX/wy2;
327	}
328
329	void AliH2F::SetSmoothRange(Float_t dx,Float_t dy)
330	{
331	Float_t wx = (fXaxis.GetXmax()-fXaxis.GetXmin()) / Float_t(fXaxis.GetNbins());
332	Float_t wy = (fYaxis.GetXmax()-fYaxis.GetXmin()) / Float_t(fYaxis.GetNbins());
333	fdx = Int_t(dx/wx);
334	fdy = Int_t(dy/wy);
335	}
336
337	TMatrix * AliH2F::GetFitMatrix(Int_t irow, Int_t icolumn)
338	{
339	SmoothCell(irow,icolumn);
340	return &fFitMatrix;
341	}
342
343	void AliH2F::SmoothCell(Int_t irow, Int_t icolumn)
344	{
345	//calculate interpolation around point irow icollumn
346	//fitting by surface of second oreder
347	Int_t index;
348	Float_t x,y,z,x2,x3,x4,y2,y3,y4;
349	Float_t sz = 0.0,szx = 0.0 ,szy = 0.0 ,szx2 = 0.0,szy2 = 0.0;
350	Float_t sx = 0.,sx2 = 0.,sx3 = 0.,sx4 = 0.;
351	Float_t sy = 0.,sy2 = 0. ,sy3 = 0. ,sy4 = 0.;
352	Float_t sxy = 0.,sxy2 = 0. ,sx2y =0.0 ,sx2y2 =0.0;
353	Float_t w;
354	Float_t sumweight = 0;
355	for (Int_t i = -fdx; i <=fdx; i++)
356	for (Int_t j = -fdy; j <=fdy; j++)
357	{
358	index = GetBin(irow+i,icolumn+j);
359	w = 1/(Float_t(ii)+fSigmaX2)1/(Float_t(j*j)+fSigmaY2);
360	z = GetBinContent(index);
361	x = irow+i;
362	x2 = xx; x3 = x2x; x4 = x2*x2;
363	y = icolumn+j;
364	y2 = yy; y3 = y2y; y4 = y2*y2;
365	sz += zw; sx += xw; sy += y*w;
366	szx += zxw; szy += zyw; szx2+= zx2w; szy2 += zy3w;
367	sx2 += x2w; sx3 += x3w; sx4 += x4*w;
368	sy2 += y2w; sy3 += y3w; sy4 += y4*w;
369	sxy += xyw; sxy2+= xy2w; sx2y+=x2yw; sx2y2+= x2y2w;
370	sumweight +=w;
371	}
372	TMatrix mat(5,5);
373	if (!fFitMatrix.IsValid()) fFitMatrix.ResizeTo(5,1);
374
375	fFitMatrix(0,0) = sz;
376	fFitMatrix(1,0) = szx;
377	fFitMatrix(2,0) = szy;
378	fFitMatrix(3,0) = szx2;
379	fFitMatrix(4,0) = szy2;
380
381	mat(0,0) = sumweight;
382	mat(0,1) = sx; mat(0,2) = sy; mat(0,3) = sx2; mat(0,4) =sy2;
383	mat(1,0) = sx; mat(1,1) = sx2; mat(1,2) = sxy; mat(1,3) = sx3; mat(1,4) =sxy2;
384	mat(2,0) = sy; mat(2,1) = sxy; mat(2,2) = sy2; mat(2,3) = sx2y; mat(2,4) =sy3 ;
385	mat(3,0) = sx2; mat(3,1) = sx3; mat(3,2) = sx2y;mat(3,3) = sx4 ; mat(3,4) =sx2y2;
386	mat(4,0) = sy2; mat(4,1) = sxy2; mat(4,2) = sy3; mat(4,3) = sx2y2;mat(4,4) =sy4;
387	}