Little changes to make g++ version 3.2 compile the hough library.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughClusterTransformer.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 "AliL3HoughClusterTransformer.h"
10 #include "AliL3MemHandler.h"
11 #include "AliL3SpacePointData.h"
12 #include "AliL3Transform.h"
13 #include "AliL3DigitData.h"
14 #include "AliL3Histogram.h"
15 #include "AliL3ClustFinderNew.h"
16
17 #if GCCVERSION == 3
18 using namespace std;
19 #endif
20
21 //_____________________________________________________________
22 // AliL3HoughClusterTransformer
23 //
24 // Hough transformation class. 
25 //
26 // This class performs the Hough transform on _clusters_, and not raw data.
27 // Therefore, it first finds the clusters using the HLT cluster finder, and
28 // then uses these as input for the transform.
29
30 ClassImp(AliL3HoughClusterTransformer)
31
32 AliL3HoughClusterTransformer::AliL3HoughClusterTransformer()
33 {
34   //Default constructor
35   fParamSpace = 0;
36   fNClusters = 0;
37   fClusters = 0;
38   fMemHandler = 0;
39 #ifdef do_mc
40   fTrackID = 0;
41 #endif
42 }
43
44 AliL3HoughClusterTransformer::AliL3HoughClusterTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
45 {
46   //Normal constructor
47   fParamSpace = 0;
48   fNClusters=0;
49   fClusters=0;
50   fMemHandler=0;
51 #ifdef do_mc
52   fTrackID=0;
53 #endif
54 }
55
56 AliL3HoughClusterTransformer::~AliL3HoughClusterTransformer()
57 {
58   DeleteHistograms();
59   if(fMemHandler)
60     delete fMemHandler;
61 #ifdef do_mc
62   if(fTrackID)
63     {
64       for(Int_t i=0; i<GetNEtaSegments(); i++)
65         {
66           if(!fTrackID[i]) continue;
67           delete fTrackID[i];
68         }
69       delete [] fTrackID;
70     }
71 #endif
72 }
73
74 //void AliL3HoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
75
76 void AliL3HoughClusterTransformer::DeleteHistograms()
77 {
78   if(fParamSpace)
79     {
80       for(Int_t i=0; i<GetNEtaSegments(); i++)
81         {
82           if(!fParamSpace[i]) continue;
83           delete fParamSpace[i];
84         }
85       delete [] fParamSpace;
86     }
87 }
88
89 void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
90                                              Int_t nybin,Double_t phimin,Double_t phimax)
91 {
92   //Create the histograms (parameter space).
93   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
94   //The arguments give the range and binning; 
95   //nxbin = #bins in kappa
96   //nybin = #bins in phi0
97   //pt_min = mimium Pt of track (corresponding to maximum kappa)
98   //phi_min = mimimum phi0 (degrees)
99   //phi_max = maximum phi0 (degrees)
100     
101   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
102   Double_t torad = AliL3Transform::Pi()/180;
103   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
104 }
105
106 void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
107                                                     Int_t nybin,Double_t ymin,Double_t ymax)
108 {
109   Char_t histname[256];
110   
111   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
112   for(Int_t i=0; i<GetNEtaSegments(); i++)
113     {
114       sprintf(histname,"paramspace_%d",i);
115       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
116     }
117 #ifdef do_mc
118   Int_t ncells = (nxbin+2)*(nybin+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 #endif
124   
125 }
126
127 void AliL3HoughClusterTransformer::Reset()
128 {
129   //Reset all the histograms
130
131   if(!fParamSpace)
132     {
133       LOG(AliL3Log::kWarning,"AliL3HoughClusterTransformer::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
141   if(fMemHandler)
142     delete fMemHandler;
143   fNClusters=0;
144   fClusters=0;
145 #ifdef do_mc
146   AliL3Histogram *hist = fParamSpace[0];
147   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
148   for(Int_t i=0; i<GetNEtaSegments(); i++)
149     memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
150 #endif
151 }
152
153 Int_t AliL3HoughClusterTransformer::GetEtaIndex(Double_t eta)
154 {
155   //Return the histogram index of the corresponding eta. 
156   
157   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
158   Double_t index = (eta-GetEtaMin())/etaslice;
159   return (Int_t)index;
160 }
161
162 inline AliL3Histogram *AliL3HoughClusterTransformer::GetHistogram(Int_t eta_index)
163 {
164   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
165     return 0;
166   if(!fParamSpace[eta_index])
167     return 0;
168   return fParamSpace[eta_index];
169 }
170
171 Double_t AliL3HoughClusterTransformer::GetEta(Int_t eta_index,Int_t slice)
172 {
173   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
174   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
175   if(slice>17) eta*=-1;
176   return eta;
177 }
178
179 void AliL3HoughClusterTransformer::FindClusters()
180 {
181   //Find the clusters
182   if(!GetDataPointer())
183     {
184       cerr<<"AliL3HoughClusterTransformer::FindClusters : Zero data pointer"<<endl;
185       return;
186     }
187
188   const Int_t maxpoints=100000;
189   const Int_t pointsize = maxpoints * sizeof(AliL3SpacePointData);
190
191   fMemHandler = new AliL3MemHandler();
192   fClusters = (AliL3SpacePointData*)fMemHandler->Allocate(pointsize);
193   AliL3ClustFinderNew *cf = new AliL3ClustFinderNew();
194   cf->InitSlice(0,GetPatch(),AliL3Transform::GetFirstRow(GetPatch()),AliL3Transform::GetLastRow(GetPatch()),maxpoints);
195   cf->SetDeconv(kFALSE);
196   cf->SetOutputArray(fClusters);
197   cf->Read(1,GetDataPointer());
198   cf->ProcessDigits();
199   fNClusters = cf->GetNumberOfClusters();
200   delete cf;
201 }
202
203 void AliL3HoughClusterTransformer::TransformCircle()
204 {
205   //Transform the input data with a circle HT.
206   //The function loops over all the data, and transforms each cluster with the equations:
207   // 
208   //kappa = 2/R*sin(phi - phi0)
209   //
210   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
211   //
212   //Each cluster then transforms into a curve in the (kappa,phi0)-space. In order to find
213   //which histogram in which the pixel should be transformed, the eta-value is calcluated
214   //and the proper histogram index is found by GetEtaIndex(eta).
215
216   FindClusters();
217   if(!fClusters)
218     {
219       LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircle","Data")
220         <<"No input data "<<ENDLOG;
221       return;
222     }
223   
224   Float_t xyz[3];
225   for(Int_t i=0; i<fNClusters; i++)
226     {
227       xyz[0] = fClusters[i].fX;
228       xyz[1] = fClusters[i].fY;
229       xyz[2] = fClusters[i].fZ;
230       Double_t eta = AliL3Transform::GetEta(xyz);
231       Int_t eta_index = GetEtaIndex(eta);
232       if(eta_index < 0 || eta_index >= GetNEtaSegments())
233         continue;
234       
235       AliL3Histogram *hist = fParamSpace[eta_index];
236       
237       //Do the transformation:
238       Double_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
239       Double_t phi = AliL3Transform::GetPhi(xyz);
240       
241       //Fill the histogram along the phirange
242       for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
243         {
244           Double_t phi0 = hist->GetBinCenterY(b);
245           Double_t kappa = 2*sin(phi - phi0)/R;
246           hist->Fill(kappa,phi0,1);
247 #ifdef do_mc
248           Int_t bin = hist->FindBin(kappa,phi0);
249           for(Int_t t=0; t<3; t++)
250             {
251               Int_t label = fClusters[i].fTrackID[t];
252               if(label < 0) break;
253               UInt_t c;
254               for(c=0; c<MaxTrack; c++)
255                 if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
256                   break;
257               if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
258               fTrackID[eta_index][bin].fLabel[c] = label;
259               fTrackID[eta_index][bin].fNHits[c]++;
260             }
261 #endif
262         }
263     }
264 }
265
266 void AliL3HoughClusterTransformer::TransformCircleC(Int_t row_range)
267 {
268   //Circle transform, using combinations of every 2 points lying
269   //on different padrows and within the same etaslice.
270   
271   FindClusters();
272   if(!fClusters)
273     LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircleC","Data")
274       <<"No input data "<<ENDLOG;
275   
276   Float_t xyz[3];
277   Double_t eta,r1,phi1,r2,phi2,kappa,psi;
278   Int_t index1,index2;
279   for(Int_t i=0; i<fNClusters; i++)
280     {
281       xyz[0] = fClusters[i].fX;
282       xyz[1] = fClusters[i].fY;
283       xyz[2] = fClusters[i].fZ;
284       r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
285       phi1 = atan2(xyz[1],xyz[0]);
286       eta = AliL3Transform::GetEta(xyz);
287       index1 = GetEtaIndex(eta);
288       if(index1 < 0 || index1 >= GetNEtaSegments()) continue;
289       
290       for(Int_t j=i+1; j<fNClusters; j++)
291         {
292           if(fClusters[j].fPadRow == fClusters[i].fPadRow) continue;
293           xyz[0] = fClusters[j].fX;
294           xyz[1] = fClusters[j].fY;
295           xyz[2] = fClusters[j].fZ;
296           eta = AliL3Transform::GetEta(xyz);
297           index2 = GetEtaIndex(eta);
298           if(index1 != index2) continue;
299           
300           AliL3Histogram *hist = fParamSpace[index1];
301           
302           r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
303           phi2 = atan2(xyz[1],xyz[0]);
304           
305           psi = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
306           kappa = 2*sin(phi1-psi) / r1;
307           hist->Fill(kappa,psi,1);
308 #ifdef do_mc
309           Int_t bin = hist->FindBin(kappa,psi);
310           for(Int_t l=0; l<3; l++)
311             {
312               for(Int_t m=0; m<3; m++)
313                 {
314                   if(fClusters[i].fTrackID[l] == fClusters[j].fTrackID[m])
315                     {
316                       Int_t label = fClusters[i].fTrackID[l];
317                       if(label < 0) continue;
318                       UInt_t c;
319                       for(c=0; c<MaxTrack; c++)
320                         if(fTrackID[index1][bin].fLabel[c] == label || fTrackID[index1][bin].fNHits[c] == 0)
321                           break;
322                       if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
323                       fTrackID[index1][bin].fLabel[c] = label;
324                       fTrackID[index1][bin].fNHits[c]++;
325                     }
326                 }
327             }
328 #endif
329         }
330     }
331 }
332
333
334 Int_t AliL3HoughClusterTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
335 {
336
337 #ifdef do_mc
338   if(eta_index < 0 || eta_index > GetNEtaSegments())
339     {
340       cerr<<"AliL3HoughClusterTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
341       return -1;
342     }
343   AliL3Histogram *hist = fParamSpace[eta_index];
344   Int_t bin = hist->FindBin(kappa,psi);
345   Int_t label=-1;
346   Int_t max=0;
347   for(UInt_t i=0; i<MaxTrack; i++)
348     {
349       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
350       if(nhits == 0) break;
351       if(nhits > max)
352         {
353           max = nhits;
354           label = fTrackID[eta_index][bin].fLabel[i];
355         }
356     }
357   return label;
358 #endif
359   cout<<"AliL3HoughClusterTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
360   return -1;
361 }
362