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