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