b59a819fe5ed4305b6e740c4e8b2067de8ddd878
[u/mrichter/AliRoot.git] / TPC / AliH2F.cxx
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.3  2000/06/30 12:07:49  kowal2
19 Updated from the TPC-PreRelease branch
20
21 Revision 1.2.4.2  2000/06/26 07:39:42  kowal2
22 Changes to obey the coding rules
23
24 Revision 1.2.4.1  2000/06/25 08:38:41  kowal2
25 Splitted from AliTPCtracking
26
27 Revision 1.2  2000/04/17 09:37:33  kowal2
28 removed obsolete AliTPCDigitsDisplay.C
29
30 Revision 1.1.4.2  2000/04/10 11:32:37  kowal2
31
32 "ROOT"-based class with some extra functionality
33
34 */
35 //
36
37 //----------------------------------------------------------------------------
38 //  Author:   Marian Ivanov
39 //
40 //  Implementation of class AliH2F
41 //
42 //-----------------------------------------------------------------------------
43
44 #include "AliH2F.h"
45 #include "TClonesArray.h"
46 #include "AliTPC.h"
47 #include "TRandom.h"
48 #include "AliTPCClusterFinder.h"
49
50
51 ClassImp(AliH2F)
52 //***********************************************************************
53 //***********************************************************************
54 //***********************************************************************
55 //***********************************************************************
56 AliH2F::AliH2F():TH2F() 
57 {
58   //
59  
60 }
61 AliH2F::AliH2F(const Text_t *name,const Text_t *title,
62                        Int_t nbinsx,Axis_t xlow,Axis_t xup
63                        ,Int_t nbinsy,Axis_t ylow,Axis_t yup):
64   TH2F(name,title,nbinsx,xlow,xup
65        ,nbinsy,ylow,yup)
66 {
67   //
68   
69 }
70      
71 AliH2F::~AliH2F() 
72 {
73   //
74 }
75
76 AliH2F::AliH2F(const AliH2F &his) 
77 {
78   //
79   
80 }
81
82 AliH2F & AliH2F::operator = (const AliH2F & his) 
83 {
84   //
85   return *this;
86 }
87
88 TClonesArray * AliH2F::FindPeaks(Float_t threshold, Float_t noise)
89 {
90   //find peaks and write it in form of AliTPCcluster to array
91     
92   //firstly we need to create object for cluster finding
93   //and fill it with contents of histogram
94   AliTPCClusterFinder cfinder;
95   cfinder.SetThreshold(threshold);
96   cfinder.SetNoise(noise);
97   cfinder.GetHisto(this);
98   return cfinder.FindPeaks3();
99 }
100
101 void AliH2F::ClearSpectrum()
102 {
103   //clera histogram
104   Int_t dimx =  fXaxis.GetNbins();
105   Int_t dimy =  fYaxis.GetNbins();
106   for (Int_t i = 0 ;i<dimx;i++)
107     for (Int_t j = 0 ;j<dimy;j++) 
108       {
109         SetCellContent(i,j,0);
110         SetCellError(i,j,0);
111       }
112 }
113
114
115 void AliH2F::AddNoise(Float_t sn)
116 {
117   // add gauss noise with sigma sn
118   Int_t dimx =  fXaxis.GetNbins();
119   Int_t dimy =  fYaxis.GetNbins();
120   for (Int_t i = 0 ;i<dimx;i++)
121     for (Int_t j = 0 ;j<dimy;j++) 
122       {
123         Float_t noise = gRandom->Gaus(0,sn);
124         Float_t oldv  =GetCellContent(i,j);
125         Float_t olds  =GetCellError(i,j);
126         if (noise >0)
127           {
128             SetCellContent(i,j,noise+oldv);
129             SetCellError(i,j,TMath::Sqrt((noise*noise+olds*olds)));
130           }
131       }
132 }
133 void AliH2F::AddGauss(Float_t x, Float_t y, 
134                           Float_t sx, Float_t sy, Float_t max)
135 {  
136   //transform to histogram coordinata  
137   Int_t dimx =  fXaxis.GetNbins();
138   Int_t dimy =  fYaxis.GetNbins();
139   Float_t dx =(GetXaxis()->GetXmax()-GetXaxis()->GetXmin())/Float_t(dimx);
140   Float_t dy =(GetYaxis()->GetXmax()-GetYaxis()->GetXmin())/Float_t(dimy);  
141   //  x=(x-GetXaxis()->GetXmin())/dx;
142   //y=(y-GetYaxis()->GetXmin())/dy;
143   sx/=dx;
144   sy/=dy;
145
146   
147   for (Int_t i = 0 ;i<dimx;i++)
148     for (Int_t j = 0 ;j<dimy;j++) 
149       {
150         Float_t x2 =GetXaxis()->GetBinCenter(i+1);
151         Float_t y2 =GetYaxis()->GetBinCenter(j+1);
152         Float_t dx2 = (x2-x)*(x2-x);
153         Float_t dy2 = (y2-y)*(y2-y);
154         Float_t amp =max*exp(-(dx2/(2*sx*sx)+dy2/(2*sy*sy)));
155         //Float_t oldv  =GetCellContent(i+1,j+1);
156         //      SetCellContent(i+1,j+1,amp+oldv);
157         Fill(x2,y2,amp);
158       }
159 }
160
161 void AliH2F::ClearUnderTh(Int_t threshold)
162 {
163   //clear histogram for bin under threshold
164   Int_t dimx =  fXaxis.GetNbins();
165   Int_t dimy =  fYaxis.GetNbins();
166   for (Int_t i = 0 ;i<=dimx;i++)
167     for (Int_t j = 0 ;j<=dimy;j++) 
168       { 
169         Float_t oldv  =GetCellContent(i,j);
170         if (oldv <threshold)
171           SetCellContent(i,j,0);
172       }
173 }
174
175 void AliH2F::Round()
176 {
177   //round float to integer 
178   Int_t dimx =  fXaxis.GetNbins();
179   Int_t dimy =  fYaxis.GetNbins();
180   for (Int_t i = 0 ;i<=dimx;i++)
181     for (Int_t j = 0 ;j<=dimy;j++) 
182       { 
183         Float_t oldv  =GetCellContent(i,j);
184         oldv=(Int_t)oldv;
185         SetCellContent(i,j,oldv);
186       }
187 }
188
189
190
191 AliH2F *AliH2F::GetSubrange2d(Float_t xmin, Float_t xmax, 
192                                       Float_t ymin, Float_t ymax)
193 {
194   //this function return pointer to the new created 
195   //histogram which is subhistogram of the 
196   //calculate number
197   //subhistogram range must be inside histogram
198
199   if (xmax<=xmin) {
200     xmin=fXaxis.GetXmin();
201     xmax=fXaxis.GetXmax();
202   }
203   if (ymax<=ymin) {
204      ymin=fYaxis.GetXmin();
205      ymax=fYaxis.GetXmax();
206   }
207
208   Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  * 
209                    Float_t(fXaxis.GetNbins()));
210   Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  * 
211                    Float_t(fYaxis.GetNbins()));
212   TString  t1 = fName ;
213   TString  t2 = fTitle ;
214   t1+="_subrange";
215   t2+="_subrange";
216   const Text_t *ktt1 = t1;
217   const Text_t *ktt2 = t2;
218   
219   AliH2F * sub = new AliH2F(ktt1,ktt2,nx,xmin,xmax,ny,ymin,ymax); 
220   
221   Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
222                     (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
223   Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
224                     (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
225   for (Int_t i=0;i<nx;i++)
226     for (Int_t j=0;j<ny;j++)
227       {
228         Int_t index1 = GetBin(i1+i,i2+j);
229         //        Int_t index2 = sub->GetBin(i,j);
230         Float_t val = GetBinContent(index1);
231         //        sub->SetBinContent(index2,val);
232         //        Float_t err = GetBinError(index1);
233         //sub->SetBinError(index2,GetBinError(index1));
234         sub->SetCellContent(i,j,val);
235       }  
236    return sub;
237 }
238
239 TH1F *AliH2F::GetAmplitudes(Float_t zmin, Float_t zmax, Float_t th, Float_t xmin, Float_t xmax, 
240                                       Float_t ymin, Float_t ymax)
241 {
242   //this function return pointer to the new created 
243   //histogram which is subhistogram of the 
244   //calculate number
245   //subhistogram range must be inside histogram
246  
247   if (xmax<=xmin) {
248     xmin=fXaxis.GetXmin();
249     xmax=fXaxis.GetXmax();
250   }
251   if (ymax<=ymin) {
252      ymin=fYaxis.GetXmin();
253      ymax=fYaxis.GetXmax();
254   }
255   Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  * 
256                    Float_t(fXaxis.GetNbins()));
257   Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  * 
258                    Float_t(fYaxis.GetNbins()));
259   TString  t1 = fName ;
260   TString  t2 = fTitle ;
261   t1+="_amplitudes";
262   t2+="_amplitudes";
263   const  Text_t *ktt1 = t1;
264   const Text_t *ktt2 = t2;
265   
266   TH1F * h = new TH1F(ktt1,ktt2,100,zmin,zmax); 
267   
268   Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
269                     (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
270   Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
271                     (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
272   for (Int_t i=0;i<nx;i++)
273     for (Int_t j=0;j<ny;j++)
274       {
275         Int_t index1 = GetBin(i1+i,i2+j);
276         Float_t val = GetBinContent(index1);
277         if (val>th) h->Fill(val);
278       }  
279    return h;
280 }
281
282 Float_t   AliH2F::GetOccupancy(Float_t th , Float_t xmin, Float_t xmax, 
283                              Float_t ymin, Float_t ymax)
284 {
285   //this function return pointer to the new created 
286   //histogram which is subhistogram of the 
287   //calculate number
288   //subhistogram range must be inside histogram
289  
290   if (xmax<=xmin) {
291     xmin=fXaxis.GetXmin();
292     xmax=fXaxis.GetXmax();
293   }
294   if (ymax<=ymin) {
295      ymin=fYaxis.GetXmin();
296      ymax=fYaxis.GetXmax();
297   }
298   Int_t nx = Int_t((xmax-xmin)/(fXaxis.GetXmax()-fXaxis.GetXmin())  * 
299                    Float_t(fXaxis.GetNbins()));
300   Int_t ny = Int_t((ymax-ymin)/(fYaxis.GetXmax()-fYaxis.GetXmin())  * 
301                    Float_t(fYaxis.GetNbins()));
302  
303   Int_t over =0; 
304   Int_t i1 = Int_t( Float_t(fXaxis.GetNbins())*(xmin-fXaxis.GetXmin())/
305                     (fXaxis.GetXmax()-fXaxis.GetXmin()) ) ;
306   Int_t i2 = Int_t( Float_t(fYaxis.GetNbins())*(ymin-fYaxis.GetXmin())/
307                     (fYaxis.GetXmax()-fYaxis.GetXmin()) ) ;
308   for (Int_t i=0;i<nx;i++)
309     for (Int_t j=0;j<ny;j++)
310       {
311         Int_t index1 = GetBin(i1+i,i2+j);
312         Float_t val = GetBinContent(index1);
313         if (val>th) over++;
314       }  
315   Int_t  all = nx*ny;
316   if (all>0)  return Float_t(over)/Float_t(all);
317   else 
318     return 0;
319 }