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