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