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