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