Added new function AliL3Hough::MergeEtaSlices which merges tracks which
[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;
44c7f8de 31 fDoMC = kFALSE;;
1eef4efc 32#ifdef do_mc
33 fTrackID = 0;
34#endif
4de874d1 35}
36
44c7f8de 37AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoMC) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
52a2a604 38{
3ef466c5 39 //Normal constructor
237d3f5c 40 fParamSpace = 0;
44c7f8de 41 if(DoMC)
42 {
43 if(patch==0)
44 fDoMC = kTRUE;
45 else
46 fDoMC = kFALSE;
47 }
1eef4efc 48#ifdef do_mc
49 fTrackID = 0;
50#endif
52a2a604 51}
52
4de874d1 53AliL3HoughTransformer::~AliL3HoughTransformer()
54{
99e7186b 55 DeleteHistograms();
1eef4efc 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
4cafa5fc 67}
68
7646f3c3 69//void AliL3HoughTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
70
99e7186b 71void AliL3HoughTransformer::DeleteHistograms()
4cafa5fc 72{
99e7186b 73 if(!fParamSpace)
74 return;
c52cf5d8 75 for(Int_t i=0; i<GetNEtaSegments(); i++)
4fc9a6a4 76 {
77 if(!fParamSpace[i]) continue;
78 delete fParamSpace[i];
79 }
99e7186b 80 delete [] fParamSpace;
4cafa5fc 81}
82
e26acabd 83void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
84 Int_t nybin,Double_t phimin,Double_t phimax)
85{
3ef466c5 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
08a66445 95 Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
26abc209 96 Double_t torad = AliL3Transform::Pi()/180;
97 CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
e26acabd 98}
99
99e7186b 100void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
101 Int_t nybin,Double_t ymin,Double_t ymax)
4cafa5fc 102{
4fc9a6a4 103
c52cf5d8 104 fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
4cafa5fc 105
99e7186b 106 Char_t histname[256];
c52cf5d8 107 for(Int_t i=0; i<GetNEtaSegments(); i++)
4de874d1 108 {
99e7186b 109 sprintf(histname,"paramspace_%d",i);
44c7f8de 110 //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.1,2,0.02,nybin,ymin,ymax);
99e7186b 111 fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
4de874d1 112 }
1eef4efc 113
114#ifdef do_mc
c980a70c 115 if(fDoMC)
116 {
88905d4d 117 AliL3Histogram *hist = fParamSpace[0];
118 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
c980a70c 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 }
1eef4efc 124#endif
4cafa5fc 125}
126
e26acabd 127void AliL3HoughTransformer::Reset()
128{
48f3c46f 129 //Reset all the histograms
e26acabd 130
131 if(!fParamSpace)
132 {
133 LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
134 <<"No histograms to reset"<<ENDLOG;
135 return;
136 }
137
c52cf5d8 138 for(Int_t i=0; i<GetNEtaSegments(); i++)
e26acabd 139 fParamSpace[i]->Reset();
1eef4efc 140#ifdef do_mc
c980a70c 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 }
1eef4efc 148#endif
e26acabd 149}
150
b5a207b4 151Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
152{
3ef466c5 153 //Return the histogram index of the corresponding eta.
154
c52cf5d8 155 Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
156 Double_t index = (eta-GetEtaMin())/etaslice;
b5a207b4 157 return (Int_t)index;
158}
159
7646f3c3 160inline 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
afd8fed4 169Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice)
7646f3c3 170{
171 Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
172 Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
7646f3c3 173 return eta;
174}
175
4fc9a6a4 176void AliL3HoughTransformer::TransformCircle()
4de874d1 177{
99e7186b 178 //Transform the input data with a circle HT.
3ef466c5 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
52a2a604 189
c52cf5d8 190 AliL3DigitRowData *tempPt = GetDataPointer();
191 if(!tempPt)
4de874d1 192 {
237d3f5c 193 LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
194 <<"No input data "<<ENDLOG;
195 return;
196 }
99e7186b 197
3ef466c5 198 //Loop over the padrows:
26abc209 199 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
4cafa5fc 200 {
3ef466c5 201 //Get the data on this padrow:
99e7186b 202 AliL3DigitData *digPt = tempPt->fDigitData;
203 if(i != (Int_t)tempPt->fRow)
4cafa5fc 204 {
d96f6a4a 205 cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
99e7186b 206 continue;
207 }
44c7f8de 208
3ef466c5 209 //Loop over the data on this padrow:
99e7186b 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;
3bb06991 215 if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
99e7186b 216 continue;
44c7f8de 217
4cafa5fc 218 Int_t sector,row;
99e7186b 219 Float_t xyz[3];
3ef466c5 220
221 //Transform data to local cartesian coordinates:
4ab9f8f0 222 AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
223 AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
3ef466c5 224
225 //Calculate the eta:
4ab9f8f0 226 Double_t eta = AliL3Transform::GetEta(xyz);
3ef466c5 227
228 //Get the corresponding index, which determines which histogram to fill:
c52cf5d8 229 Int_t eta_index = GetEtaIndex(eta);
230 if(eta_index < 0 || eta_index >= GetNEtaSegments())
99e7186b 231 continue;
4cafa5fc 232
e26acabd 233 //Get the correct histogrampointer:
99e7186b 234 AliL3Histogram *hist = fParamSpace[eta_index];
235 if(!hist)
4cafa5fc 236 {
4fc9a6a4 237 printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
99e7186b 238 continue;
4cafa5fc 239 }
4cafa5fc 240
3ef466c5 241 //Do the transformation:
242 Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
4ab9f8f0 243 Float_t phi = AliL3Transform::GetPhi(xyz);
4de874d1 244
99e7186b 245 //Fill the histogram along the phirange
246 for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
4de874d1 247 {
99e7186b 248 Float_t phi0 = hist->GetBinCenterY(b);
249 Float_t kappa = 2*sin(phi - phi0)/R;
44c7f8de 250 hist->Fill(kappa,phi0,1);//charge);
1eef4efc 251#ifdef do_mc
c980a70c 252 if(fDoMC)
1eef4efc 253 {
c980a70c 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 }
1eef4efc 267 }
268#endif
4de874d1 269 }
4de874d1 270 }
3ef466c5 271
272 //Move the data pointer to the next padrow:
48f3c46f 273 AliL3MemHandler::UpdateRowPointer(tempPt);
4de874d1 274 }
4de874d1 275}
276
3ef466c5 277void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
b1886074 278{
279 //Circle transform, using combinations of every 2 points lying
280 //on different padrows and within the same etaslice.
281
c52cf5d8 282 AliL3DigitRowData *tempPt = GetDataPointer();
b1886074 283 if(!tempPt)
237d3f5c 284 LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
285 <<"No input data "<<ENDLOG;
286
3ef466c5 287 Int_t counter=0;
26abc209 288 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
b1886074 289 {
3ef466c5 290 counter += tempPt->fNDigit;
b1886074 291 AliL3MemHandler::UpdateRowPointer(tempPt);
292 }
3ef466c5 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;
1eef4efc 300 Int_t trackID[3];
3ef466c5 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;
b1886074 308 Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
b1886074 309 Float_t xyz[3];
3ef466c5 310
311 counter=0;
c52cf5d8 312 tempPt = GetDataPointer();
3ef466c5 313
26abc209 314 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
b1886074 315 {
3ef466c5 316 AliL3DigitData *digPt = tempPt->fDigitData;
317 for(UInt_t di=0; di<tempPt->fNDigit; di++)
b1886074 318 {
3ef466c5 319 charge = digPt[di].fCharge;
b1886074 320 pad = digPt[di].fPad;
321 time = digPt[di].fTime;
4ab9f8f0 322 AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
323 AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
324 eta = AliL3Transform::GetEta(xyz);
3ef466c5 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;
93657878 330#ifdef do_mc
c980a70c 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 }
93657878 337#endif
3ef466c5 338 counter++;
339 }
340 AliL3MemHandler::UpdateRowPointer(tempPt);
341 }
342
343 for(Int_t i=0; i<total_digits; i++)
344 {
c52cf5d8 345 if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
3ef466c5 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;
b1886074 353
354 //Get the correct histogrampointer:
3ef466c5 355 AliL3Histogram *hist = fParamSpace[ind];
b1886074 356 if(!hist)
357 {
3ef466c5 358 printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
b1886074 359 continue;
360 }
3ef466c5 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;
44c7f8de 369 hist->Fill(kappa,phi_0,1);//tot_charge);
1eef4efc 370#ifdef do_mc
c980a70c 371 if(fDoMC)
1eef4efc 372 {
c980a70c 373 Int_t bin = hist->FindBin(kappa,phi_0);
374 for(Int_t l=0; l<3; l++)
1eef4efc 375 {
c980a70c 376 for(Int_t m=0; m<3; m++)
1eef4efc 377 {
c980a70c 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 }
1eef4efc 390 }
391 }
392 }
393#endif
b1886074 394 }
395 }
3ef466c5 396 delete [] digits;
b1886074 397}
398
4fc9a6a4 399void AliL3HoughTransformer::TransformLine()
400{
401 //Do a line transform on the data.
402
403
c52cf5d8 404 AliL3DigitRowData *tempPt = GetDataPointer();
405 if(!tempPt)
4fc9a6a4 406 {
237d3f5c 407 LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
408 <<"No input data "<<ENDLOG;
409 return;
410 }
26abc209 411
412 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
4fc9a6a4 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;
3bb06991 425 if(charge < GetLowerThreshold())
4fc9a6a4 426 continue;
427 Int_t sector,row;
428 Float_t xyz[3];
4ab9f8f0 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);
b5a207b4 432 Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
c52cf5d8 433 if(eta_index < 0 || eta_index >= GetNEtaSegments())
4fc9a6a4 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 }
48f3c46f 450 AliL3MemHandler::UpdateRowPointer(tempPt);
4fc9a6a4 451 }
452
453}
1eef4efc 454
455Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
456{
c980a70c 457 if(!fDoMC)
458 {
459 cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
460 return -1;
461 }
462
1eef4efc 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