Little changes to make it work with new base class.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformer.cxx
1 //$Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ASV 
5
6 #include "AliL3MemHandler.h"
7 #include "AliL3Logging.h"
8 #include "AliL3HoughTransformer.h"
9 #include "AliL3Transform.h"
10 #include "AliL3DigitData.h"
11 #include "AliL3Histogram.h"
12
13 //_____________________________________________________________
14 // AliL3HoughTransformer
15 //
16 // Hough transformation class
17 //
18
19 ClassImp(AliL3HoughTransformer)
20
21 AliL3HoughTransformer::AliL3HoughTransformer()
22 {
23   //Default constructor
24   fParamSpace = 0;
25 }
26
27 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
28 {
29   //Normal constructor
30   fParamSpace = 0;
31 }
32
33 AliL3HoughTransformer::~AliL3HoughTransformer()
34 {
35   DeleteHistograms();
36 }
37
38 //void AliL3HoughTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
39
40 void AliL3HoughTransformer::DeleteHistograms()
41 {
42   if(!fParamSpace)
43     return;
44   for(Int_t i=0; i<GetNEtaSegments(); i++)
45     {
46       if(!fParamSpace[i]) continue;
47       delete fParamSpace[i];
48     }
49   delete [] fParamSpace;
50 }
51
52 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
53                                              Int_t nybin,Double_t phimin,Double_t phimax)
54 {
55   //Create the histograms (parameter space).
56   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
57   //The arguments give the range and binning; 
58   //nxbin = #bins in kappa
59   //nybin = #bins in phi0
60   //pt_min = mimium Pt of track (corresponding to maximum kappa)
61   //phi_min = mimimum phi0 (degrees)
62   //phi_max = maximum phi0 (degrees)
63     
64   Double_t bfact = 0.0029980;
65   Double_t bfield = 0.2;
66   Double_t x = bfact*bfield/pt_min;
67   Double_t torad = AliL3Transform::Pi()/180;
68   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
69 }
70
71 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
72                                              Int_t nybin,Double_t ymin,Double_t ymax)
73 {
74   
75   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
76   
77   Char_t histname[256];
78   for(Int_t i=0; i<GetNEtaSegments(); i++)
79     {
80       sprintf(histname,"paramspace_%d",i);
81       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
82     }
83 }
84
85 void AliL3HoughTransformer::Reset()
86 {
87   //Reset all the histograms
88
89   if(!fParamSpace)
90     {
91       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
92         <<"No histograms to reset"<<ENDLOG;
93       return;
94     }
95   
96   for(Int_t i=0; i<GetNEtaSegments(); i++)
97     fParamSpace[i]->Reset();
98 }
99
100 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
101 {
102   //Return the histogram index of the corresponding eta. 
103
104   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
105   Double_t index = (eta-GetEtaMin())/etaslice;
106   return (Int_t)index;
107 }
108
109 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
110 {
111   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
112     return 0;
113   if(!fParamSpace[eta_index])
114     return 0;
115   return fParamSpace[eta_index];
116 }
117
118 Double_t AliL3HoughTransformer::GetEta(Int_t eta_index)
119 {
120   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
121   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
122   if(GetSlice()>17) eta*=-1;
123   return eta;
124 }
125
126 void AliL3HoughTransformer::TransformCircle()
127 {
128   //Transform the input data with a circle HT.
129   //The function loops over all the data, and transforms each pixel with the equations:
130   // 
131   //kappa = 2/R*sin(phi - phi0)
132   //
133   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
134   //
135   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
136   //which histogram in which the pixel should be transformed, the eta-value is calcluated
137   //and the proper histogram index is found by GetEtaIndex(eta).
138
139
140   AliL3DigitRowData *tempPt = GetDataPointer();
141   if(!tempPt)
142     {
143       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
144         <<"No input data "<<ENDLOG;
145       return;
146     }
147   
148   //Loop over the padrows:
149   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
150     {
151       //Get the data on this padrow:
152       AliL3DigitData *digPt = tempPt->fDigitData;
153       if(i != (Int_t)tempPt->fRow)
154         {
155           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
156           continue;
157         }
158
159       //Loop over the data on this padrow:
160       for(UInt_t j=0; j<tempPt->fNDigit; j++)
161         {
162           UShort_t charge = digPt[j].fCharge;
163           UChar_t pad = digPt[j].fPad;
164           UShort_t time = digPt[j].fTime;
165           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
166             continue;
167           Int_t sector,row;
168           Float_t xyz[3];
169           
170           //Transform data to local cartesian coordinates:
171           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
172           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
173           
174           //Calculate the eta:
175           Double_t eta = AliL3Transform::GetEta(xyz);
176           
177           //Get the corresponding index, which determines which histogram to fill:
178           Int_t eta_index = GetEtaIndex(eta);
179           if(eta_index < 0 || eta_index >= GetNEtaSegments())
180             continue;
181           
182           //Get the correct histogrampointer:
183           AliL3Histogram *hist = fParamSpace[eta_index];
184           if(!hist)
185             {
186               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
187               continue;
188             }
189
190           //Do the transformation:
191           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
192           Float_t phi = AliL3Transform::GetPhi(xyz);
193           
194           //Fill the histogram along the phirange
195           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
196             {
197               Float_t phi0 = hist->GetBinCenterY(b);
198               Float_t kappa = 2*sin(phi - phi0)/R;
199               hist->Fill(kappa,phi0,charge);
200             }
201         }
202       
203       //Move the data pointer to the next padrow:
204       AliL3MemHandler::UpdateRowPointer(tempPt);
205     }
206 }
207
208 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
209 {
210   //Circle transform, using combinations of every 2 points lying
211   //on different padrows and within the same etaslice.
212   
213   AliL3DigitRowData *tempPt = GetDataPointer();
214   if(!tempPt)
215     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
216       <<"No input data "<<ENDLOG;
217  
218   Int_t counter=0;
219   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
220     {
221       counter += tempPt->fNDigit;
222       AliL3MemHandler::UpdateRowPointer(tempPt);
223     }
224   
225   struct Digit {
226     Int_t row;
227     Double_t r;
228     Double_t phi;
229     Int_t eta_index;
230     Int_t charge;
231   };
232   
233   Digit *digits = new Digit[counter];
234   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
235   
236   Int_t total_digits=counter;
237   Int_t sector,row,tot_charge,pad,time,charge;
238   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
239   Float_t xyz[3];
240   
241   counter=0;
242   tempPt = GetDataPointer();
243   
244   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
245     {
246       AliL3DigitData *digPt = tempPt->fDigitData;
247       for(UInt_t di=0; di<tempPt->fNDigit; di++)
248         {
249           charge = digPt[di].fCharge;
250           pad = digPt[di].fPad;
251           time = digPt[di].fTime;
252           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
253           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
254           eta = AliL3Transform::GetEta(xyz);
255           digits[counter].row = i;
256           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
257           digits[counter].phi = atan2(xyz[1],xyz[0]);
258           digits[counter].eta_index = GetEtaIndex(eta);
259           digits[counter].charge = charge;
260           counter++;
261         }
262       AliL3MemHandler::UpdateRowPointer(tempPt);
263     }
264   
265   for(Int_t i=0; i<total_digits; i++)
266     {
267       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
268       Int_t ind = digits[i].eta_index;
269       
270       for(Int_t j=i+1; j<total_digits; j++)
271         {
272           if(digits[i].row == digits[j].row) continue;
273           if(digits[i].eta_index != digits[j].eta_index) continue;
274           if(digits[i].row + row_range < digits[j].row) break;
275           
276           //Get the correct histogrampointer:
277           AliL3Histogram *hist = fParamSpace[ind];
278           if(!hist)
279             {
280               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
281               continue;
282             }
283           
284           r1 = digits[i].r;
285           phi1 = digits[i].phi;
286           r2 = digits[j].r;
287           phi2 = digits[j].phi;
288           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
289           kappa = 2*sin(phi2-phi_0)/r2;
290           tot_charge = digits[i].charge + digits[j].charge;
291           hist->Fill(kappa,phi_0,tot_charge);
292         }
293     }
294   delete [] digits;
295 }
296
297 void AliL3HoughTransformer::TransformLine()
298 {
299   //Do a line transform on the data.
300
301   
302   AliL3DigitRowData *tempPt = GetDataPointer();
303   if(!tempPt)
304     {
305       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
306         <<"No input data "<<ENDLOG;
307       return;
308     }
309     
310   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
311     {
312       AliL3DigitData *digPt = tempPt->fDigitData;
313       if(i != (Int_t)tempPt->fRow)
314         {
315           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
316           continue;
317         }
318       for(UInt_t j=0; j<tempPt->fNDigit; j++)
319         {
320           UShort_t charge = digPt[j].fCharge;
321           UChar_t pad = digPt[j].fPad;
322           UShort_t time = digPt[j].fTime;
323           if(charge < GetLowerThreshold())
324             continue;
325           Int_t sector,row;
326           Float_t xyz[3];
327           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
328           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
329           Float_t eta = AliL3Transform::GetEta(xyz);
330           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
331           if(eta_index < 0 || eta_index >= GetNEtaSegments())
332             continue;
333           
334           //Get the correct histogram:
335           AliL3Histogram *hist = fParamSpace[eta_index];
336           if(!hist)
337             {
338               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
339               continue;
340             }
341           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
342             {
343               Double_t theta = hist->GetBinCenterX(xbin);
344               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
345               hist->Fill(theta,rho,charge);
346             }
347         }
348       AliL3MemHandler::UpdateRowPointer(tempPt);
349     }
350   
351 }