aebd5a1e645208f022676f5bca3b901b9c72acbd
[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 "AliL3HistogramAdaptive.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   fDoMC = kFALSE;
26 #ifdef do_mc
27   fTrackID = 0;
28 #endif
29 }
30
31 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
32 {
33   //Normal constructor
34   fParamSpace = 0;
35   fDoMC = kFALSE;
36 #ifdef do_mc
37   fTrackID = 0;
38 #endif
39 }
40
41 AliL3HoughTransformer::~AliL3HoughTransformer()
42 {
43   DeleteHistograms();
44 #ifdef do_mc
45   if(fTrackID)
46     {
47       for(Int_t i=0; i<GetNEtaSegments(); i++)
48         {
49           if(!fTrackID[i]) continue;
50           delete fTrackID[i];
51         }
52       delete [] fTrackID;
53     }
54 #endif
55 }
56
57 //void AliL3HoughTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
58
59 void AliL3HoughTransformer::DeleteHistograms()
60 {
61   if(!fParamSpace)
62     return;
63   for(Int_t i=0; i<GetNEtaSegments(); i++)
64     {
65       if(!fParamSpace[i]) continue;
66       delete fParamSpace[i];
67     }
68   delete [] fParamSpace;
69 }
70
71 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
72                                              Int_t nybin,Double_t phimin,Double_t phimax)
73 {
74   //Create the histograms (parameter space).
75   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
76   //The arguments give the range and binning; 
77   //nxbin = #bins in kappa
78   //nybin = #bins in phi0
79   //pt_min = mimium Pt of track (corresponding to maximum kappa)
80   //phi_min = mimimum phi0 (degrees)
81   //phi_max = maximum phi0 (degrees)
82     
83   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
84   Double_t torad = AliL3Transform::Pi()/180;
85   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
86 }
87
88 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
89                                              Int_t nybin,Double_t ymin,Double_t ymax)
90 {
91   
92   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
93   
94   Char_t histname[256];
95   for(Int_t i=0; i<GetNEtaSegments(); i++)
96     {
97       sprintf(histname,"paramspace_%d",i);
98       //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.1,1,0.05,64,ymin,ymax);
99       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
100     }
101   
102 #ifdef do_mc
103   if(fDoMC)
104     {
105       AliL3Histogram *hist = fParamSpace[0];
106       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
107       cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
108       fTrackID = new TrackIndex*[GetNEtaSegments()];
109       for(Int_t i=0; i<GetNEtaSegments(); i++)
110         fTrackID[i] = new TrackIndex[ncells];
111     }
112 #endif
113 }
114
115 void AliL3HoughTransformer::Reset()
116 {
117   //Reset all the histograms
118
119   if(!fParamSpace)
120     {
121       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
122         <<"No histograms to reset"<<ENDLOG;
123       return;
124     }
125   
126   for(Int_t i=0; i<GetNEtaSegments(); i++)
127     fParamSpace[i]->Reset();
128 #ifdef do_mc
129   if(fDoMC)
130     {
131       AliL3Histogram *hist = fParamSpace[0];
132       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
133       for(Int_t i=0; i<GetNEtaSegments(); i++)
134         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
135     }
136 #endif
137 }
138
139 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
140 {
141   //Return the histogram index of the corresponding eta. 
142
143   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
144   Double_t index = (eta-GetEtaMin())/etaslice;
145   return (Int_t)index;
146 }
147
148 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
149 {
150   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
151     return 0;
152   if(!fParamSpace[eta_index])
153     return 0;
154   return fParamSpace[eta_index];
155 }
156
157 Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice)
158 {
159   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
160   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
161   if(slice>17) eta*=-1;
162   return eta;
163 }
164
165 void AliL3HoughTransformer::TransformCircle()
166 {
167   //Transform the input data with a circle HT.
168   //The function loops over all the data, and transforms each pixel with the equations:
169   // 
170   //kappa = 2/R*sin(phi - phi0)
171   //
172   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
173   //
174   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
175   //which histogram in which the pixel should be transformed, the eta-value is calcluated
176   //and the proper histogram index is found by GetEtaIndex(eta).
177
178
179   AliL3DigitRowData *tempPt = GetDataPointer();
180   if(!tempPt)
181     {
182       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
183         <<"No input data "<<ENDLOG;
184       return;
185     }
186   
187   //Loop over the padrows:
188   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
189     {
190       //Get the data on this padrow:
191       AliL3DigitData *digPt = tempPt->fDigitData;
192       if(i != (Int_t)tempPt->fRow)
193         {
194           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
195           continue;
196         }
197
198       //Loop over the data on this padrow:
199       for(UInt_t j=0; j<tempPt->fNDigit; j++)
200         {
201           UShort_t charge = digPt[j].fCharge;
202           UChar_t pad = digPt[j].fPad;
203           UShort_t time = digPt[j].fTime;
204           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
205             continue;
206           Int_t sector,row;
207           Float_t xyz[3];
208           
209           //Transform data to local cartesian coordinates:
210           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
211           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
212           
213           //Calculate the eta:
214           Double_t eta = AliL3Transform::GetEta(xyz);
215           
216           //Get the corresponding index, which determines which histogram to fill:
217           Int_t eta_index = GetEtaIndex(eta);
218           if(eta_index < 0 || eta_index >= GetNEtaSegments())
219             continue;
220           
221           //Get the correct histogrampointer:
222           AliL3Histogram *hist = fParamSpace[eta_index];
223           if(!hist)
224             {
225               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
226               continue;
227             }
228
229           //Do the transformation:
230           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
231           Float_t phi = AliL3Transform::GetPhi(xyz);
232           
233           //Fill the histogram along the phirange
234           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
235             {
236               Float_t phi0 = hist->GetBinCenterY(b);
237               Float_t kappa = 2*sin(phi - phi0)/R;
238               hist->Fill(kappa,phi0,charge);
239 #ifdef do_mc
240               if(fDoMC)
241                 {
242                   Int_t bin = hist->FindBin(kappa,phi0);
243                   for(Int_t t=0; t<3; t++)
244                     {
245                       Int_t label = digPt[j].fTrackID[t];
246                       if(label < 0) break;
247                       UInt_t c;
248                       for(c=0; c<MaxTrack; c++)
249                         if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
250                           break;
251                       if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
252                       fTrackID[eta_index][bin].fLabel[c] = label;
253                       fTrackID[eta_index][bin].fNHits[c]++;
254                     }
255                 }
256 #endif
257             }
258         }
259       
260       //Move the data pointer to the next padrow:
261       AliL3MemHandler::UpdateRowPointer(tempPt);
262     }
263 }
264
265 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
266 {
267   //Circle transform, using combinations of every 2 points lying
268   //on different padrows and within the same etaslice.
269   
270   AliL3DigitRowData *tempPt = GetDataPointer();
271   if(!tempPt)
272     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
273       <<"No input data "<<ENDLOG;
274  
275   Int_t counter=0;
276   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
277     {
278       counter += tempPt->fNDigit;
279       AliL3MemHandler::UpdateRowPointer(tempPt);
280     }
281   
282   struct Digit {
283     Int_t row;
284     Double_t r;
285     Double_t phi;
286     Int_t eta_index;
287     Int_t charge;
288     Int_t trackID[3];
289   };
290   
291   Digit *digits = new Digit[counter];
292   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
293   
294   Int_t total_digits=counter;
295   Int_t sector,row,tot_charge,pad,time,charge;
296   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
297   Float_t xyz[3];
298   
299   counter=0;
300   tempPt = GetDataPointer();
301   
302   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
303     {
304       AliL3DigitData *digPt = tempPt->fDigitData;
305       for(UInt_t di=0; di<tempPt->fNDigit; di++)
306         {
307           charge = digPt[di].fCharge;
308           pad = digPt[di].fPad;
309           time = digPt[di].fTime;
310           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
311           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
312           eta = AliL3Transform::GetEta(xyz);
313           digits[counter].row = i;
314           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
315           digits[counter].phi = atan2(xyz[1],xyz[0]);
316           digits[counter].eta_index = GetEtaIndex(eta);
317           digits[counter].charge = charge;
318 #ifdef do_mc
319           if(fDoMC)
320             {
321               digits[counter].trackID[0] = digPt[di].fTrackID[0];
322               digits[counter].trackID[1] = digPt[di].fTrackID[1];
323               digits[counter].trackID[2] = digPt[di].fTrackID[2];
324             }
325 #endif
326           counter++;
327         }
328       AliL3MemHandler::UpdateRowPointer(tempPt);
329     }
330   
331   for(Int_t i=0; i<total_digits; i++)
332     {
333       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
334       Int_t ind = digits[i].eta_index;
335       
336       for(Int_t j=i+1; j<total_digits; j++)
337         {
338           if(digits[i].row == digits[j].row) continue;
339           if(digits[i].eta_index != digits[j].eta_index) continue;
340           if(digits[i].row + row_range < digits[j].row) break;
341           
342           //Get the correct histogrampointer:
343           AliL3Histogram *hist = fParamSpace[ind];
344           if(!hist)
345             {
346               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
347               continue;
348             }
349           
350           r1 = digits[i].r;
351           phi1 = digits[i].phi;
352           r2 = digits[j].r;
353           phi2 = digits[j].phi;
354           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
355           kappa = 2*sin(phi2-phi_0)/r2;
356           tot_charge = digits[i].charge + digits[j].charge;
357           hist->Fill(kappa,phi_0,tot_charge);
358 #ifdef do_mc
359           if(fDoMC)
360             {
361               Int_t bin = hist->FindBin(kappa,phi_0);
362               for(Int_t l=0; l<3; l++)
363                 {
364                   for(Int_t m=0; m<3; m++)
365                     {
366                       if(digits[i].trackID[l] == digits[j].trackID[m])
367                         {
368                           Int_t label = digits[i].trackID[l];
369                           if(label < 0) continue;
370                           UInt_t c;
371                           for(c=0; c<MaxTrack; c++)
372                             if(fTrackID[ind][bin].fLabel[c] == label || fTrackID[ind][bin].fNHits[c] == 0)
373                               break;
374                           if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
375                           fTrackID[ind][bin].fLabel[c] = label;
376                           fTrackID[ind][bin].fNHits[c]++;
377                         }
378                     }
379                 }
380             }
381 #endif
382         }
383     }
384   delete [] digits;
385 }
386
387 void AliL3HoughTransformer::TransformLine()
388 {
389   //Do a line transform on the data.
390
391   
392   AliL3DigitRowData *tempPt = GetDataPointer();
393   if(!tempPt)
394     {
395       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
396         <<"No input data "<<ENDLOG;
397       return;
398     }
399     
400   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
401     {
402       AliL3DigitData *digPt = tempPt->fDigitData;
403       if(i != (Int_t)tempPt->fRow)
404         {
405           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
406           continue;
407         }
408       for(UInt_t j=0; j<tempPt->fNDigit; j++)
409         {
410           UShort_t charge = digPt[j].fCharge;
411           UChar_t pad = digPt[j].fPad;
412           UShort_t time = digPt[j].fTime;
413           if(charge < GetLowerThreshold())
414             continue;
415           Int_t sector,row;
416           Float_t xyz[3];
417           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
418           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
419           Float_t eta = AliL3Transform::GetEta(xyz);
420           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
421           if(eta_index < 0 || eta_index >= GetNEtaSegments())
422             continue;
423           
424           //Get the correct histogram:
425           AliL3Histogram *hist = fParamSpace[eta_index];
426           if(!hist)
427             {
428               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
429               continue;
430             }
431           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
432             {
433               Double_t theta = hist->GetBinCenterX(xbin);
434               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
435               hist->Fill(theta,rho,charge);
436             }
437         }
438       AliL3MemHandler::UpdateRowPointer(tempPt);
439     }
440   
441 }
442
443 Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
444 {
445   if(!fDoMC)
446     {
447       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
448       return -1;
449     }
450   
451 #ifdef do_mc
452   if(eta_index < 0 || eta_index > GetNEtaSegments())
453     {
454       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
455       return -1;
456     }
457   AliL3Histogram *hist = fParamSpace[eta_index];
458   Int_t bin = hist->FindBin(kappa,psi);
459   Int_t label=-1;
460   Int_t max=0;
461   for(UInt_t i=0; i<MaxTrack; i++)
462     {
463       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
464       if(nhits == 0) break;
465       if(nhits > max)
466         {
467           max = nhits;
468           label = fTrackID[eta_index][bin].fLabel[i];
469         }
470     }
471   return label;
472 #endif
473   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
474   return -1;
475 }
476