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