Made a new abstract base class; AliL3HoughBaseTransformer for different implementations
[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   
32    fParamSpace = 0;
33 }
34
35 AliL3HoughTransformer::~AliL3HoughTransformer()
36 {
37   DeleteHistograms();
38 }
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   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*ToRad,phimax*ToRad);
68 }
69
70 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
71                                              Int_t nybin,Double_t ymin,Double_t ymax)
72 {
73   
74   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
75   
76   Char_t histname[256];
77   for(Int_t i=0; i<GetNEtaSegments(); i++)
78     {
79       sprintf(histname,"paramspace_%d",i);
80       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
81     }
82 }
83
84 void AliL3HoughTransformer::Reset()
85 {
86   //Reset all the histograms
87
88   if(!fParamSpace)
89     {
90       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
91         <<"No histograms to reset"<<ENDLOG;
92       return;
93     }
94   
95   for(Int_t i=0; i<GetNEtaSegments(); i++)
96     fParamSpace[i]->Reset();
97 }
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 void AliL3HoughTransformer::TransformCircle()
110 {
111   //Transform the input data with a circle HT.
112   //The function loops over all the data, and transforms each pixel with the equations:
113   // 
114   //kappa = 2/R*sin(phi - phi0)
115   //
116   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
117   //
118   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
119   //which histogram in which the pixel should be transformed, the eta-value is calcluated
120   //and the proper histogram index is found by GetEtaIndex(eta).
121
122
123   AliL3DigitRowData *tempPt = GetDataPointer();
124   if(!tempPt)
125     {
126       printf("\nAliL3HoughTransformer::TransformCircle : No input data!!!\n\n");
127       return;
128     }
129   
130   //Loop over the padrows:
131   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
132     {
133       //Get the data on this padrow:
134       AliL3DigitData *digPt = tempPt->fDigitData;
135       if(i != (Int_t)tempPt->fRow)
136         {
137           printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
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(charge <= GetThreshold())
148             continue;
149           Int_t sector,row;
150           Float_t xyz[3];
151           
152           //Transform data to local cartesian coordinates:
153           fTransform->Slice2Sector(GetSlice(),i,sector,row);
154           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
155           
156           //Calculate the eta:
157           Double_t eta = fTransform->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 = fTransform->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     printf("\nAliL3HoughTransformer::TransformCircleC() : Zero data pointer\n");
198   
199   Int_t counter=0;
200   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
201     {
202       counter += tempPt->fNDigit;
203       AliL3MemHandler::UpdateRowPointer(tempPt);
204     }
205   
206   struct Digit {
207     Int_t row;
208     Double_t r;
209     Double_t phi;
210     Int_t eta_index;
211     Int_t charge;
212   };
213   
214   Digit *digits = new Digit[counter];
215   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
216   
217   Int_t total_digits=counter;
218   Int_t sector,row,tot_charge,pad,time,charge;
219   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
220   Float_t xyz[3];
221   
222   counter=0;
223   tempPt = GetDataPointer();
224   
225   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
226     {
227       AliL3DigitData *digPt = tempPt->fDigitData;
228       for(UInt_t di=0; di<tempPt->fNDigit; di++)
229         {
230           charge = digPt[di].fCharge;
231           pad = digPt[di].fPad;
232           time = digPt[di].fTime;
233           fTransform->Slice2Sector(GetSlice(),i,sector,row);
234           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
235           eta = fTransform->GetEta(xyz);
236           digits[counter].row = i;
237           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
238           digits[counter].phi = atan2(xyz[1],xyz[0]);
239           digits[counter].eta_index = GetEtaIndex(eta);
240           digits[counter].charge = charge;
241           counter++;
242         }
243       AliL3MemHandler::UpdateRowPointer(tempPt);
244     }
245   
246   for(Int_t i=0; i<total_digits; i++)
247     {
248       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
249       Int_t ind = digits[i].eta_index;
250       
251       for(Int_t j=i+1; j<total_digits; j++)
252         {
253           if(digits[i].row == digits[j].row) continue;
254           if(digits[i].eta_index != digits[j].eta_index) continue;
255           if(digits[i].row + row_range < digits[j].row) break;
256           
257           //Get the correct histogrampointer:
258           AliL3Histogram *hist = fParamSpace[ind];
259           if(!hist)
260             {
261               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
262               continue;
263             }
264           
265           r1 = digits[i].r;
266           phi1 = digits[i].phi;
267           r2 = digits[j].r;
268           phi2 = digits[j].phi;
269           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
270           kappa = 2*sin(phi2-phi_0)/r2;
271           tot_charge = digits[i].charge + digits[j].charge;
272           hist->Fill(kappa,phi_0,tot_charge);
273         }
274     }
275   delete [] digits;
276 }
277
278 void AliL3HoughTransformer::TransformLine()
279 {
280   //Do a line transform on the data.
281
282   
283   AliL3DigitRowData *tempPt = GetDataPointer();
284   if(!tempPt)
285     {
286       printf("\nAliL3HoughTransformer::TransformLine : No input data!!!\n\n");
287       return;
288     }
289   
290   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
291     {
292       AliL3DigitData *digPt = tempPt->fDigitData;
293       if(i != (Int_t)tempPt->fRow)
294         {
295           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
296           continue;
297         }
298       for(UInt_t j=0; j<tempPt->fNDigit; j++)
299         {
300           UShort_t charge = digPt[j].fCharge;
301           UChar_t pad = digPt[j].fPad;
302           UShort_t time = digPt[j].fTime;
303           if(charge < GetThreshold())
304             continue;
305           Int_t sector,row;
306           Float_t xyz[3];
307           fTransform->Slice2Sector(GetSlice(),i,sector,row);
308           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
309           Float_t eta = fTransform->GetEta(xyz);
310           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
311           if(eta_index < 0 || eta_index >= GetNEtaSegments())
312             continue;
313           
314           //Get the correct histogram:
315           AliL3Histogram *hist = fParamSpace[eta_index];
316           if(!hist)
317             {
318               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
319               continue;
320             }
321           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
322             {
323               Double_t theta = hist->GetBinCenterX(xbin);
324               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
325               hist->Fill(theta,rho,charge);
326             }
327         }
328       AliL3MemHandler::UpdateRowPointer(tempPt);
329     }
330   
331 }