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