24e427a7c5f4e9530cfdfa07af370125e86f91cd
[u/mrichter/AliRoot.git] / HLT / TPCLib / tracking / AliHLTTPCHistogramAdaptive.cxx
1 // @(#) $Id: AliHLTTPCHistogramAdaptive.cxx,v 1.1 2006/11/30 17:45:43 hristov Exp 
2 // origin: hough/AliL3HistogramAdaptive.cxx,v 1.13 Thu Jun 23 17:46:54 2005 UTC by hristov
3
4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
5 //*-- Copyright &copy ALICE HLT Group
6
7 #include "AliHLTStdIncludes.h"
8 #include "AliHLTTPCLogging.h"
9 #include "AliHLTTPCHistogramAdaptive.h"
10 #include "AliHLTTPCTransform.h"
11 #include "AliHLTTPCTrack.h"
12
13 #if __GNUC__ >= 3
14 using namespace std;
15 #endif
16
17 //_____________________________________________________________
18 // AliHLTTPCHistogramAdaptive
19 //
20 // 2D histogram class adapted for kappa and psi as used in the Circle Hough Transform.
21 // The bins in kappa is not linear, but has a width which is specified by argument
22 // ptres in the constructor. This gives the relative pt resolution which should
23 // be kept throughout the kappa range. 
24
25 ClassImp(AliHLTTPCHistogramAdaptive)
26
27 AliHLTTPCHistogramAdaptive::AliHLTTPCHistogramAdaptive() : AliHLTTPCHistogram()
28 {
29   //default ctor
30   fKappaBins=0;
31 }
32
33   
34 AliHLTTPCHistogramAdaptive::AliHLTTPCHistogramAdaptive(Char_t *name,Double_t minpt,Double_t maxpt,Double_t ptres,
35                                                Int_t nybins,Double_t ymin,Double_t ymax)
36 {
37   //normal ctor
38   strcpy(fName,name);
39   
40   fPtres = ptres;
41   fXmin = -1*AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/minpt;
42   fXmax = AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/minpt;
43
44   fMinPt = minpt;
45   fMaxPt = maxpt;
46   fNxbins = InitKappaBins();
47   fNybins = nybins;
48   
49   fYmin = ymin;
50   fYmax = ymax;
51   fFirstXbin=1;
52   fFirstYbin=1;
53   fLastXbin = fNxbins;
54   fLastYbin = fNybins;
55   fNcells = (fNxbins+2)*(fNybins+2);
56
57   fThreshold=0;
58   fContent = new Int_t[fNcells];
59   Reset();
60 }
61
62 AliHLTTPCHistogramAdaptive::~AliHLTTPCHistogramAdaptive()
63 {
64   //dtor
65   if(fKappaBins)
66     delete [] fKappaBins;
67 }
68
69 Int_t AliHLTTPCHistogramAdaptive::InitKappaBins()
70 {
71   //Here a LUT for the kappa values created. This has to be done since
72   //the binwidth in kappa is not constant, but change according to the
73   //set relative resolution in pt.
74   //Since the kappa values are symmetric about origo, the size of the
75   //LUT is half of the total number of bins in kappa direction.
76   
77   Double_t pt = fMinPt,deltapt,localpt;
78   Int_t bin=0;
79   
80   while(pt < fMaxPt)
81     {
82       localpt = pt;
83       deltapt = fPtres*localpt*localpt;
84       pt += 2*deltapt;
85       bin++;
86     }
87   fKappaBins = new Double_t[bin+1];
88   pt=fMinPt;
89   bin=0;
90   fKappaBins[bin] = AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/fMinPt; 
91   while(pt < fMaxPt)
92     {
93       localpt = pt;
94       deltapt = fPtres*localpt*localpt;
95       pt += 2*deltapt;                      //*2 because pt +- 1/2*deltapt is one bin
96       bin++;
97       fKappaBins[bin] = AliHLTTPCTransform::GetBFact()*AliHLTTPCTransform::GetBField()/pt;
98     }
99   return (bin+1)*2; //Both negative and positive kappa.
100 }
101
102
103 void AliHLTTPCHistogramAdaptive::Fill(Double_t x,Double_t y,Int_t weight)
104 {
105   //Fill a given bin in the histogram
106   Int_t bin = FindBin(x,y);
107   if(bin < 0)
108     return;
109   AddBinContent(bin,weight);
110
111 }
112
113 Int_t AliHLTTPCHistogramAdaptive::FindBin(Double_t x,Double_t y) const
114 {
115   //Find a bin in the histogram  
116   Int_t xbin = FindXbin(x);
117   Int_t ybin = FindYbin(y);
118   
119   if(!xbin || !ybin) 
120     return -1;
121   return GetBin(xbin,ybin);
122 }
123
124 Int_t AliHLTTPCHistogramAdaptive::FindXbin(Double_t x) const
125 {
126   //Find X bin in the histogram
127   if(x < fXmin || x > fXmax || fabs(x) < fKappaBins[(fNxbins/2-1)])
128     return 0;
129   
130   //Remember that kappa value is decreasing with bin number!
131   //Also, the bin numbering starts at 1 and ends at fNxbins,
132   //so the corresponding elements in the LUT is bin - 1.
133
134   Int_t bin=0;
135   while(bin < fNxbins/2)
136     {
137       if(fabs(x) <= fKappaBins[bin] && fabs(x) > fKappaBins[bin+1])
138         break;
139       bin++;
140     }
141   if(x < 0)
142     return bin + 1;
143   else 
144     return fNxbins - bin;
145   
146 }
147
148 Int_t AliHLTTPCHistogramAdaptive::FindYbin(Double_t y) const
149 {
150   //Find Y bin in the histogram
151   if(y < fYmin || y > fYmax)
152     return 0;
153   
154   return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
155 }
156
157 Double_t AliHLTTPCHistogramAdaptive::GetBinCenterX(Int_t xbin) const
158 {
159   //Returns bin center in X
160   if(xbin < fFirstXbin || xbin > fLastXbin)
161     {
162       LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHistogramAdaptive::GetBinCenterX","Bin-value")
163         <<"XBinvalue out of range "<<xbin<<ENDLOG;
164       return 0;
165     }
166   
167   //The bin numbers go from 1 to fNxbins, so the corresponding
168   //element in the LUT is xbin - 1. This is the reason why we 
169   //substract a 1 here:
170   
171   Int_t bin = xbin;
172   bin -= 1;
173   if(bin >= fNxbins/2)
174     bin = fNxbins - 1 - bin;
175   
176   //Remember again that the kappa-values are _decreasing_ with bin number.
177   
178   Double_t binwidth = fKappaBins[bin] - fKappaBins[bin+1];
179   Double_t kappa = fKappaBins[bin] - 0.5*binwidth;
180   if(xbin < fNxbins/2)
181     return -1.*kappa;
182   else
183     return kappa;
184
185 }
186
187 Double_t AliHLTTPCHistogramAdaptive::GetBinCenterY(Int_t ybin) const
188 {
189   //Returns bin center in Y
190   if(ybin < fFirstYbin || ybin > fLastYbin)
191     {
192       LOG(AliHLTTPCLog::kError,"AliHLTTPCHistogramAdaptive::GetBinCenterY","ybin")
193         <<"Bin-value out of range "<<ybin<<ENDLOG;
194       return -1;
195     }
196   Double_t binwidth = (fYmax - fYmin) / fNybins;
197   return fYmin + (ybin-0.5) * binwidth;
198
199 }
200
201
202 void AliHLTTPCHistogramAdaptive::Draw(Char_t *option)
203 {
204   //Draw the histogram
205 #ifdef use_root
206   if(!fRootHisto)
207     CreateRootHisto();
208   
209   Double_t kappa,psi;
210   Int_t content,bin;
211   for(Int_t i=fFirstXbin; i<=fLastXbin; i++)
212     {
213       kappa = GetBinCenterX(i);
214       for(Int_t j=fFirstYbin; j<=fLastYbin; j++)
215         {
216           psi = GetBinCenterY(j);
217           bin = GetBin(i,j);
218           content = GetBinContent(bin);
219           fRootHisto->Fill(kappa,psi,content);
220         }
221     }
222   fRootHisto->Draw(option);
223   return;
224 #else
225   cerr<<"AliHLTTPCHistogramAdaptive::Draw : You need to compile with ROOT in order to draw histogram"<<endl;
226 #endif
227 }
228
229 void AliHLTTPCHistogramAdaptive::Print() const
230 {
231   //Print the contents of the histogram
232   cout<<"Printing content of histogram "<<fName<<endl;
233   for(Int_t i=0; i<fNcells; i++)
234     {
235       if(GetBinContent(i)==0) continue;
236       cout<<"Bin "<<i<<": "<<GetBinContent(i)<<endl;
237     }
238
239 }