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