removed obsolete AliTPCDigitsDisplay.C
[u/mrichter/AliRoot.git] / TPC / AliH2F.cxx
CommitLineData
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$
18Revision 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
42ClassImp(AliH2F)
43//***********************************************************************
44//***********************************************************************
45//***********************************************************************
46//***********************************************************************
47AliH2F::AliH2F():TH2F()
48{
49 fSigmaX2 = 1.;
50 fSigmaY2 = 1.;
51 fdx = 1;
52 fdy=1;
53 fFitMatrix.ResizeTo(5,1);
54}
55AliH2F::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
68AliH2F::~AliH2F()
69{
70}
71
72AliH2F::AliH2F(const AliH2F &)
73{
74}
75
76AliH2F & AliH2F::operator = (const AliH2F &)
77{
78 return *this;
79}
80
81TClonesArray * 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
94void 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
107void 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((noise*noise+olds*olds)));
121 }
122 }
123}
124void 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 =max*exp(-(dx2/(2*sx*sx)+dy2/(2*sy*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
152void 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
165void 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
180AliH2F *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
228TH1F *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
271Float_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
319void 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
329void 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
337TMatrix * AliH2F::GetFitMatrix(Int_t irow, Int_t icolumn)
338{
339 SmoothCell(irow,icolumn);
340 return &fFitMatrix;
341}
342
343void 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(i*i)+fSigmaX2)*1/(Float_t(j*j)+fSigmaY2);
360 z = GetBinContent(index);
361 x = irow+i;
362 x2 = x*x; x3 = x2*x; x4 = x2*x2;
363 y = icolumn+j;
364 y2 = y*y; y3 = y2*y; y4 = y2*y2;
365 sz += z*w; sx += x*w; sy += y*w;
366 szx += z*x*w; szy += z*y*w; szx2+= z*x2*w; szy2 += z*y3*w;
367 sx2 += x2*w; sx3 += x3*w; sx4 += x4*w;
368 sy2 += y2*w; sy3 += y3*w; sy4 += y4*w;
369 sxy += x*y*w; sxy2+= x*y2*w; sx2y+=x2*y*w; sx2y2+= x2*y2*w;
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}