Some additional changes related to the previous changes. AliL3Transform
[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   fParamSpace = 0;
26 }
27
28 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
29 {
30   //Normal constructor
31   fParamSpace = 0;
32 }
33
34 AliL3HoughTransformer::~AliL3HoughTransformer()
35 {
36   DeleteHistograms();
37 }
38
39 void AliL3HoughTransformer::DeleteHistograms()
40 {
41   if(!fParamSpace)
42     return;
43   for(Int_t i=0; i<GetNEtaSegments(); i++)
44     {
45       if(!fParamSpace[i]) continue;
46       delete fParamSpace[i];
47     }
48   delete [] fParamSpace;
49 }
50
51 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
52                                              Int_t nybin,Double_t phimin,Double_t phimax)
53 {
54   //Create the histograms (parameter space).
55   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
56   //The arguments give the range and binning; 
57   //nxbin = #bins in kappa
58   //nybin = #bins in phi0
59   //pt_min = mimium Pt of track (corresponding to maximum kappa)
60   //phi_min = mimimum phi0 (degrees)
61   //phi_max = maximum phi0 (degrees)
62     
63   Double_t bfact = 0.0029980;
64   Double_t bfield = 0.2;
65   Double_t x = bfact*bfield/pt_min;
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   if(!fTransform)
130     {
131       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Transformer")
132         <<"No AliL3Transform object"<<ENDLOG;
133       return;
134     }
135   
136   //Loop over the padrows:
137   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
138     {
139       //Get the data on this padrow:
140       AliL3DigitData *digPt = tempPt->fDigitData;
141       if(i != (Int_t)tempPt->fRow)
142         {
143           printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
144           continue;
145         }
146       
147       //Loop over the data on this padrow:
148       for(UInt_t j=0; j<tempPt->fNDigit; j++)
149         {
150           UShort_t charge = digPt[j].fCharge;
151           UChar_t pad = digPt[j].fPad;
152           UShort_t time = digPt[j].fTime;
153           if(charge <= GetThreshold())
154             continue;
155           Int_t sector,row;
156           Float_t xyz[3];
157           
158           //Transform data to local cartesian coordinates:
159           fTransform->Slice2Sector(GetSlice(),i,sector,row);
160           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
161           
162           //Calculate the eta:
163           Double_t eta = fTransform->GetEta(xyz);
164           
165           //Get the corresponding index, which determines which histogram to fill:
166           Int_t eta_index = GetEtaIndex(eta);
167           if(eta_index < 0 || eta_index >= GetNEtaSegments())
168             continue;
169           
170           //Get the correct histogrampointer:
171           AliL3Histogram *hist = fParamSpace[eta_index];
172           if(!hist)
173             {
174               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
175               continue;
176             }
177
178           //Do the transformation:
179           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
180           Float_t phi = fTransform->GetPhi(xyz);
181           
182           //Fill the histogram along the phirange
183           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
184             {
185               Float_t phi0 = hist->GetBinCenterY(b);
186               Float_t kappa = 2*sin(phi - phi0)/R;
187               hist->Fill(kappa,phi0,charge);
188             }
189         }
190       
191       //Move the data pointer to the next padrow:
192       AliL3MemHandler::UpdateRowPointer(tempPt);
193     }
194 }
195
196 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
197 {
198   //Circle transform, using combinations of every 2 points lying
199   //on different padrows and within the same etaslice.
200   
201   AliL3DigitRowData *tempPt = GetDataPointer();
202   if(!tempPt)
203     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
204       <<"No input data "<<ENDLOG;
205  
206   if(!fTransform)
207     {
208       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Transformer")
209         <<"No AliL3Transform object"<<ENDLOG;
210       return;
211     }
212   
213   Int_t counter=0;
214   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
215     {
216       counter += tempPt->fNDigit;
217       AliL3MemHandler::UpdateRowPointer(tempPt);
218     }
219   
220   struct Digit {
221     Int_t row;
222     Double_t r;
223     Double_t phi;
224     Int_t eta_index;
225     Int_t charge;
226   };
227   
228   Digit *digits = new Digit[counter];
229   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
230   
231   Int_t total_digits=counter;
232   Int_t sector,row,tot_charge,pad,time,charge;
233   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
234   Float_t xyz[3];
235   
236   counter=0;
237   tempPt = GetDataPointer();
238   
239   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
240     {
241       AliL3DigitData *digPt = tempPt->fDigitData;
242       for(UInt_t di=0; di<tempPt->fNDigit; di++)
243         {
244           charge = digPt[di].fCharge;
245           pad = digPt[di].fPad;
246           time = digPt[di].fTime;
247           fTransform->Slice2Sector(GetSlice(),i,sector,row);
248           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
249           eta = fTransform->GetEta(xyz);
250           digits[counter].row = i;
251           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
252           digits[counter].phi = atan2(xyz[1],xyz[0]);
253           digits[counter].eta_index = GetEtaIndex(eta);
254           digits[counter].charge = charge;
255           counter++;
256         }
257       AliL3MemHandler::UpdateRowPointer(tempPt);
258     }
259   
260   for(Int_t i=0; i<total_digits; i++)
261     {
262       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
263       Int_t ind = digits[i].eta_index;
264       
265       for(Int_t j=i+1; j<total_digits; j++)
266         {
267           if(digits[i].row == digits[j].row) continue;
268           if(digits[i].eta_index != digits[j].eta_index) continue;
269           if(digits[i].row + row_range < digits[j].row) break;
270           
271           //Get the correct histogrampointer:
272           AliL3Histogram *hist = fParamSpace[ind];
273           if(!hist)
274             {
275               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
276               continue;
277             }
278           
279           r1 = digits[i].r;
280           phi1 = digits[i].phi;
281           r2 = digits[j].r;
282           phi2 = digits[j].phi;
283           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
284           kappa = 2*sin(phi2-phi_0)/r2;
285           tot_charge = digits[i].charge + digits[j].charge;
286           hist->Fill(kappa,phi_0,tot_charge);
287         }
288     }
289   delete [] digits;
290 }
291
292 void AliL3HoughTransformer::TransformLine()
293 {
294   //Do a line transform on the data.
295
296   
297   AliL3DigitRowData *tempPt = GetDataPointer();
298   if(!tempPt)
299     {
300       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
301         <<"No input data "<<ENDLOG;
302       return;
303     }
304   if(!fTransform)
305     {
306       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Transformer")
307         <<"No AliL3Transform object"<<ENDLOG;
308       return;
309     }
310   
311   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
312     {
313       AliL3DigitData *digPt = tempPt->fDigitData;
314       if(i != (Int_t)tempPt->fRow)
315         {
316           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
317           continue;
318         }
319       for(UInt_t j=0; j<tempPt->fNDigit; j++)
320         {
321           UShort_t charge = digPt[j].fCharge;
322           UChar_t pad = digPt[j].fPad;
323           UShort_t time = digPt[j].fTime;
324           if(charge < GetThreshold())
325             continue;
326           Int_t sector,row;
327           Float_t xyz[3];
328           fTransform->Slice2Sector(GetSlice(),i,sector,row);
329           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
330           Float_t eta = fTransform->GetEta(xyz);
331           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
332           if(eta_index < 0 || eta_index >= GetNEtaSegments())
333             continue;
334           
335           //Get the correct histogram:
336           AliL3Histogram *hist = fParamSpace[eta_index];
337           if(!hist)
338             {
339               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
340               continue;
341             }
342           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
343             {
344               Double_t theta = hist->GetBinCenterX(xbin);
345               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
346               hist->Fill(theta,rho,charge);
347             }
348         }
349       AliL3MemHandler::UpdateRowPointer(tempPt);
350     }
351   
352 }