--- /dev/null
+/**************************************************************************
+ * 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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff