Fixed missing initialization for some variables at the edges of the hough space
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerRow.cxx
CommitLineData
0bd0c1ef 1// @(#) $Id$
2
de3c3890 3
0bd0c1ef 4// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
5//*-- Copyright &copy ALICE HLT Group
6
7#include "AliL3StandardIncludes.h"
8
9#include "AliL3Logging.h"
10#include "AliL3MemHandler.h"
11#include "AliL3Transform.h"
12#include "AliL3DigitData.h"
13#include "AliL3HistogramAdaptive.h"
de3c3890 14#include "AliL3HoughTrack.h"
0bd0c1ef 15#include "AliL3HoughTransformerRow.h"
16
de3c3890 17#if GCCVERSION == 3
0bd0c1ef 18using namespace std;
19#endif
20
21//_____________________________________________________________
22// AliL3HoughTransformerRow
23//
24// TPC rows Hough transformation class
25//
26
27ClassImp(AliL3HoughTransformerRow)
28
917e711b 29UChar_t **AliL3HoughTransformerRow::fgRowCount = 0;
30UChar_t **AliL3HoughTransformerRow::fgGapCount = 0;
31UChar_t **AliL3HoughTransformerRow::fgCurrentRowCount = 0;
0bd0c1ef 32#ifdef do_mc
bd2f8772 33AliL3TrackIndex **AliL3HoughTransformerRow::fgTrackID = 0;
0bd0c1ef 34#endif
917e711b 35UChar_t *AliL3HoughTransformerRow::fgTrackNRows = 0;
36UChar_t *AliL3HoughTransformerRow::fgTrackFirstRow = 0;
37UChar_t *AliL3HoughTransformerRow::fgTrackLastRow = 0;
0bd0c1ef 38
917e711b 39Float_t AliL3HoughTransformerRow::fgBeta1 = AliL3Transform::Row2X(79)/pow(AliL3Transform::Row2X(79),2);
40Float_t AliL3HoughTransformerRow::fgBeta2 = (AliL3Transform::Row2X(158)+6.0)/pow((AliL3Transform::Row2X(158)+6.0),2);
0bd0c1ef 41
42AliL3HoughTransformerRow::AliL3HoughTransformerRow()
43{
44 //Default constructor
45 fParamSpace = 0;
46 fDoMC = kFALSE;;
47
0bd0c1ef 48 fLUTforwardZ=0;
49 fLUTforwardZ2=0;
50 fLUTbackwardZ=0;
51 fLUTbackwardZ2=0;
52}
53
917e711b 54AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,netasegments,zvertex)
0bd0c1ef 55{
56 //Normal constructor
57 fParamSpace = 0;
58 fDoMC = kFALSE;
59 fDoMC=kFALSE;
60#ifdef do_mc
61 fDoMC = kTRUE;
62#endif
63
0bd0c1ef 64 fLUTforwardZ=0;
65 fLUTforwardZ2=0;
66 fLUTbackwardZ=0;
67 fLUTbackwardZ2=0;
68}
69
70AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
71{
917e711b 72 //Destructor
0bd0c1ef 73 DeleteHistograms();
74#ifdef do_mc
917e711b 75 if(fgTrackID)
0bd0c1ef 76 {
77 for(Int_t i=0; i<GetNEtaSegments(); i++)
78 {
917e711b 79 if(!fgTrackID[i]) continue;
74e77d1b 80 delete [] fgTrackID[i];
0bd0c1ef 81 }
917e711b 82 delete [] fgTrackID;
83 fgTrackID = 0;
0bd0c1ef 84 }
85#endif
86
917e711b 87 if(fgRowCount)
0bd0c1ef 88 {
89 for(Int_t i=0; i<GetNEtaSegments(); i++)
90 {
917e711b 91 if(!fgRowCount[i]) continue;
92 delete [] fgRowCount[i];
0bd0c1ef 93 }
917e711b 94 delete [] fgRowCount;
95 fgRowCount = 0;
0bd0c1ef 96 }
917e711b 97 if(fgGapCount)
0bd0c1ef 98 {
99 for(Int_t i=0; i<GetNEtaSegments(); i++)
100 {
917e711b 101 if(!fgGapCount[i]) continue;
102 delete [] fgGapCount[i];
0bd0c1ef 103 }
917e711b 104 delete [] fgGapCount;
105 fgGapCount = 0;
0bd0c1ef 106 }
917e711b 107 if(fgCurrentRowCount)
0bd0c1ef 108 {
109 for(Int_t i=0; i<GetNEtaSegments(); i++)
110 {
917e711b 111 if(fgCurrentRowCount[i])
112 delete [] fgCurrentRowCount[i];
0bd0c1ef 113 }
917e711b 114 delete [] fgCurrentRowCount;
115 fgCurrentRowCount = 0;
0bd0c1ef 116 }
917e711b 117 if(fgTrackNRows)
de3c3890 118 {
917e711b 119 delete [] fgTrackNRows;
120 fgTrackNRows = 0;
de3c3890 121 }
917e711b 122 if(fgTrackFirstRow)
de3c3890 123 {
917e711b 124 delete [] fgTrackFirstRow;
125 fgTrackFirstRow = 0;
de3c3890 126 }
917e711b 127 if(fgTrackLastRow)
de3c3890 128 {
917e711b 129 delete [] fgTrackLastRow;
130 fgTrackLastRow = 0;
de3c3890 131 }
0bd0c1ef 132}
133
134void AliL3HoughTransformerRow::DeleteHistograms()
135{
917e711b 136 // Clean up
0bd0c1ef 137 delete[] fLUTforwardZ;
138 delete[] fLUTforwardZ2;
139 delete[] fLUTbackwardZ;
140 delete[] fLUTbackwardZ2;
141
142 fLUTforwardZ=0;
143 fLUTforwardZ2=0;
144 fLUTbackwardZ=0;
145 fLUTbackwardZ2=0;
146
147 if(!fParamSpace)
148 return;
149 for(Int_t i=0; i<GetNEtaSegments(); i++)
150 {
151 if(!fParamSpace[i]) continue;
152 delete fParamSpace[i];
153 }
154 delete [] fParamSpace;
155}
156
0bd0c1ef 157void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
158 Int_t nybin,Float_t ymin,Float_t ymax)
159{
917e711b 160 //Create the histograms (parameter space)
161 //nxbin = #bins in X
162 //nybin = #bins in Y
163 //xmin xmax ymin ymax = histogram limits in X and Y
0bd0c1ef 164 fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
165
166 Char_t histname[256];
167 for(Int_t i=0; i<GetNEtaSegments(); i++)
168 {
169 sprintf(histname,"paramspace_%d",i);
170 fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
171 }
172
173#ifdef do_mc
174 if(fDoMC)
175 {
176 AliL3Histogram *hist = fParamSpace[0];
e81ffe36 177 Int_t ncellsx = (hist->GetNbinsX()+3)/2;
178 Int_t ncellsy = (hist->GetNbinsY()+3)/2;
179 Int_t ncells = ncellsx*ncellsy;
917e711b 180 if(!fgTrackID)
0bd0c1ef 181 {
de3c3890 182 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
bd2f8772 183 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fgTrackID"<<ENDLOG;
184 fgTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
0bd0c1ef 185 for(Int_t i=0; i<GetNEtaSegments(); i++)
bd2f8772 186 fgTrackID[i] = new AliL3TrackIndex[ncells];
0bd0c1ef 187 }
188 }
189#endif
190 AliL3Histogram *hist = fParamSpace[0];
191 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
917e711b 192 if(!fgRowCount)
0bd0c1ef 193 {
de3c3890 194 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
917e711b 195 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgRowCount"<<ENDLOG;
196 fgRowCount = new UChar_t*[GetNEtaSegments()];
0bd0c1ef 197 for(Int_t i=0; i<GetNEtaSegments(); i++)
917e711b 198 fgRowCount[i] = new UChar_t[ncells];
0bd0c1ef 199 }
917e711b 200 if(!fgGapCount)
0bd0c1ef 201 {
de3c3890 202 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
917e711b 203 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgGapCount"<<ENDLOG;
204 fgGapCount = new UChar_t*[GetNEtaSegments()];
0bd0c1ef 205 for(Int_t i=0; i<GetNEtaSegments(); i++)
917e711b 206 fgGapCount[i] = new UChar_t[ncells];
0bd0c1ef 207 }
917e711b 208 if(!fgCurrentRowCount)
0bd0c1ef 209 {
de3c3890 210 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
917e711b 211 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgCurrentRowCount"<<ENDLOG;
212 fgCurrentRowCount = new UChar_t*[GetNEtaSegments()];
0bd0c1ef 213 for(Int_t i=0; i<GetNEtaSegments(); i++)
917e711b 214 fgCurrentRowCount[i] = new UChar_t[ncells];
0bd0c1ef 215 }
216
917e711b 217 if(!fgTrackNRows)
de3c3890 218 {
219 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
917e711b 220 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackNRows"<<ENDLOG;
221 fgTrackNRows = new UChar_t[ncells];
de3c3890 222 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
917e711b 223 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackFirstRow"<<ENDLOG;
224 fgTrackFirstRow = new UChar_t[ncells];
de3c3890 225 LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
917e711b 226 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackLastRow"<<ENDLOG;
227 fgTrackLastRow = new UChar_t[ncells];
de3c3890 228
229 AliL3HoughTrack track;
230 Int_t xmin = hist->GetFirstXbin();
231 Int_t xmax = hist->GetLastXbin();
232 Int_t ymin = hist->GetFirstYbin();
233 Int_t ymax = hist->GetLastYbin();
234 Int_t nxbins = hist->GetNbinsX()+2;
74e77d1b 235 for(Int_t ybin=ymin-1; ybin<=(ymax+1); ybin++)
de3c3890 236 {
74e77d1b 237 for(Int_t xbin=xmin-1; xbin<=(xmax+1); xbin++)
de3c3890 238 {
239 //cvetan: we get strange warning on gcc-2.95
240 //warning: large integer implicitly truncated to unsigned type
917e711b 241 fgTrackNRows[xbin + ybin*nxbins] = 99999;
de3c3890 242 for(Int_t deltay = -1; deltay <= 1; deltay += 2) {
243 for(Int_t deltax = -1; deltax <= 1; deltax += 2) {
244
245 Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)deltax);
246 Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)deltay);
247
917e711b 248 Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
249 Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
de3c3890 250 track.SetTrackParameters(kappa,psi,1);
917e711b 251 Bool_t firstrow = kFALSE;
de3c3890 252 UInt_t maxfirstrow = 0;
253 UInt_t maxlastrow = 0;
254 UInt_t curfirstrow = 0;
255 UInt_t curlastrow = 0;
256 for(Int_t j=AliL3Transform::GetFirstRow(0); j<=AliL3Transform::GetLastRow(5); j++)
257 {
258 Float_t hit[3];
259 if(!track.GetCrossingPoint(j,hit)) continue;
260 AliL3Transform::LocHLT2Raw(hit,0,j);
261 if(hit[1]>=0 && hit[1]<AliL3Transform::GetNPads(j))
262 {
917e711b 263 if(!firstrow) {
de3c3890 264 curfirstrow = j;
917e711b 265 firstrow = kTRUE;
de3c3890 266 }
267 curlastrow = j;
268 }
269 else {
917e711b 270 if(firstrow) {
271 firstrow = kFALSE;
de3c3890 272 if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
273 maxfirstrow = curfirstrow;
274 maxlastrow = curlastrow;
275 }
276 }
277 }
278 }
279 if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
280 maxfirstrow = curfirstrow;
281 maxlastrow = curlastrow;
282 }
917e711b 283 if((maxlastrow-maxfirstrow) < fgTrackNRows[xbin + ybin*nxbins]) {
284 fgTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
285 fgTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
286 fgTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
de3c3890 287 }
288 }
289 }
917e711b 290 // cout<<" fgTrackNRows "<<xbin<<" "<<ybin<<" "<<(Int_t)fgTrackNRows[xbin + ybin*nxbins]<<" "<<endl;
de3c3890 291 }
292 }
293 }
0bd0c1ef 294
295 //create lookup table for z of the digits
296 Int_t ntimebins = AliL3Transform::GetNTimeBins();
297 fLUTforwardZ = new Float_t[ntimebins];
298 fLUTforwardZ2 = new Float_t[ntimebins];
299 fLUTbackwardZ = new Float_t[ntimebins];
300 fLUTbackwardZ2 = new Float_t[ntimebins];
301 for(Int_t i=0; i<ntimebins; i++){
302 Float_t z;
303 z=AliL3Transform::GetZFast(0,i,GetZVertex());
304 fLUTforwardZ[i]=z;
305 fLUTforwardZ2[i]=z*z;
306 z=AliL3Transform::GetZFast(18,i,GetZVertex());
307 fLUTbackwardZ[i]=z;
308 fLUTbackwardZ2[i]=z*z;
309 }
310
311}
312
313void AliL3HoughTransformerRow::Reset()
314{
917e711b 315 //Reset all the histograms. Should be done when processing new slice
0bd0c1ef 316
317 if(!fParamSpace)
318 {
319 LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
320 <<"No histograms to reset"<<ENDLOG;
321 return;
322 }
323
324 for(Int_t i=0; i<GetNEtaSegments(); i++)
325 fParamSpace[i]->Reset();
326
327#ifdef do_mc
328 if(fDoMC)
329 {
330 AliL3Histogram *hist = fParamSpace[0];
e81ffe36 331 Int_t ncellsx = (hist->GetNbinsX()+3)/2;
332 Int_t ncellsy = (hist->GetNbinsY()+3)/2;
333 Int_t ncells = ncellsx*ncellsy;
0bd0c1ef 334 for(Int_t i=0; i<GetNEtaSegments(); i++)
bd2f8772 335 memset(fgTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
0bd0c1ef 336 }
337#endif
338 AliL3Histogram *hist = fParamSpace[0];
339 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
340 for(Int_t i=0; i<GetNEtaSegments(); i++)
341 {
917e711b 342 memset(fgRowCount[i],0,ncells*sizeof(UChar_t));
343 memset(fgGapCount[i],1,ncells*sizeof(UChar_t));
344 // memset(fgCurrentRowCount[i],160,ncells*sizeof(UChar_t));
345 memcpy(fgCurrentRowCount[i],fgTrackFirstRow,ncells*sizeof(UChar_t));
0bd0c1ef 346 }
347}
348
bd2f8772 349Int_t AliL3HoughTransformerRow::GetEtaIndex(Double_t eta) const
0bd0c1ef 350{
351 //Return the histogram index of the corresponding eta.
352
353 Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
354 Double_t index = (eta-GetEtaMin())/etaslice;
355 return (Int_t)index;
356}
357
917e711b 358inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t etaindex)
0bd0c1ef 359{
917e711b 360 // Return a pointer to the histogram which contains etaindex eta slice
361 if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
0bd0c1ef 362 return 0;
917e711b 363 if(!fParamSpace[etaindex])
0bd0c1ef 364 return 0;
917e711b 365 return fParamSpace[etaindex];
0bd0c1ef 366}
367
974fb714 368Double_t AliL3HoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/) const
0bd0c1ef 369{
917e711b 370 // Return eta calculated in the middle of the eta slice
371 Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
0bd0c1ef 372 Double_t eta=0;
917e711b 373 eta=(Double_t)((etaindex+0.5)*etaslice);
0bd0c1ef 374 return eta;
375}
376
917e711b 377struct AliL3EtaRow {
378 UChar_t fStartPad; //First pad in the cluster
379 UChar_t fEndPad; //Last pad in the cluster
380 Bool_t fIsFound; //Is the cluster already found
381 Float_t fStartY; //Y position of the first pad in the cluster
0bd0c1ef 382#ifdef do_mc
917e711b 383 Int_t fMcLabels[MaxTrack]; //Array to store mc labels inside cluster
0bd0c1ef 384#endif
385};
386
387void AliL3HoughTransformerRow::TransformCircle()
388{
917e711b 389 //Do the Hough Transform
390 Float_t beta1 = fgBeta1;
391 Float_t beta2 = fgBeta2;
392 Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
de3c3890 393
917e711b 394 Int_t netasegments = GetNEtaSegments();
395 Double_t etamin = GetEtaMin();
396 Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
0bd0c1ef 397
917e711b 398 Int_t lowerthreshold = GetLowerThreshold();
0bd0c1ef 399
400 //Assumes that all the etaslice histos are the same!
401 AliL3Histogram *h = fParamSpace[0];
917e711b 402 Float_t ymin = h->GetYmin();
de3c3890 403 //Float_t y_max = h->GetYmax();
917e711b 404 Float_t histbin = h->GetBinWidthY();
405 Int_t firstbin = h->GetFirstYbin();
406 Int_t lastbin = h->GetLastYbin();
407 Float_t xmin = h->GetXmin();
408 Float_t xmax = h->GetXmax();
409 Float_t xbin = (xmax-xmin)/h->GetNbinsX();
410 Int_t firstbinx = h->GetFirstXbin()-1;
411 Int_t lastbinx = h->GetLastXbin()+1;
0bd0c1ef 412 Int_t nbinx = h->GetNbinsX()+2;
413
917e711b 414 UChar_t lastpad;
415 Int_t lastetaindex;
416 AliL3EtaRow *etaclust = new AliL3EtaRow[netasegments];
0bd0c1ef 417
418 AliL3DigitRowData *tempPt = GetDataPointer();
419 if(!tempPt)
420 {
421 LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
422 <<"No input data "<<ENDLOG;
423 return;
424 }
425
426 Int_t ipatch = GetPatch();
917e711b 427 Double_t padpitch = AliL3Transform::GetPadPitchWidth(ipatch);
0bd0c1ef 428 Int_t islice = GetSlice();
429 Float_t *fLUTz;
430 Float_t *fLUTz2;
431 if(islice < 18) {
432 fLUTz = fLUTforwardZ;
433 fLUTz2 = fLUTforwardZ2;
434 }
435 else {
436 fLUTz = fLUTbackwardZ;
437 fLUTz2 = fLUTbackwardZ2;
438 }
439
0bd0c1ef 440 //Loop over the padrows:
441 for(UChar_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
442 {
443 Int_t npads = AliL3Transform::GetNPads((Int_t)i)-1;
444
917e711b 445 lastpad = 255;
0bd0c1ef 446 //Flush eta clusters array
917e711b 447 memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));
0bd0c1ef 448
449 Float_t x = AliL3Transform::Row2X((Int_t)i);
450 Float_t x2 = x*x;
dd7d3870 451 Float_t y=0,r2=0;
0bd0c1ef 452
0bd0c1ef 453 //Get the data on this padrow:
454 AliL3DigitData *digPt = tempPt->fDigitData;
455 if((Int_t)i != (Int_t)tempPt->fRow)
456 {
de3c3890 457 LOG(AliL3Log::kError,"AliL3HoughTransformerRow::TransformCircle","Data")
458 <<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<ENDLOG;
0bd0c1ef 459 continue;
460 }
461 // cout<<" Starting row "<<i<<endl;
462 //Loop over the data on this padrow:
463 for(UInt_t j=0; j<tempPt->fNDigit; j++)
464 {
465 UShort_t charge = digPt[j].fCharge;
917e711b 466 if((Int_t)charge <= lowerthreshold)
0bd0c1ef 467 continue;
468 UChar_t pad = digPt[j].fPad;
469 UShort_t time = digPt[j].fTime;
de3c3890 470
917e711b 471 if(pad != lastpad)
0bd0c1ef 472 {
917e711b 473 y = (pad-0.5*npads)*padpitch;
0bd0c1ef 474 Float_t y2 = y*y;
475 r2 = x2 + y2;
917e711b 476 lastetaindex = -1;
0bd0c1ef 477 }
478
479 //Transform data to local cartesian coordinates:
480 Float_t z = fLUTz[(Int_t)time];
481 Float_t z2 = fLUTz2[(Int_t)time];
482 //Calculate the eta:
483 Double_t r = sqrt(r2+z2);
484 Double_t eta = 0.5 * log((r+z)/(r-z));
485 //Get the corresponding index, which determines which histogram to fill:
917e711b 486 Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
0bd0c1ef 487
488#ifndef do_mc
917e711b 489 if(etaindex == lastetaindex) continue;
0bd0c1ef 490#endif
917e711b 491 // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
0bd0c1ef 492
917e711b 493 if(etaindex < 0 || etaindex >= netasegments)
0bd0c1ef 494 continue;
495
917e711b 496 if(!etaclust[etaindex].fIsFound)
0bd0c1ef 497 {
917e711b 498 etaclust[etaindex].fStartPad = pad;
499 etaclust[etaindex].fEndPad = pad;
500 etaclust[etaindex].fStartY = y - padpitch/2.0;
501 etaclust[etaindex].fIsFound = 1;
0bd0c1ef 502#ifdef do_mc
503 if(fDoMC)
504 {
505 for(Int_t t=0; t<3; t++)
506 {
507 Int_t label = digPt[j].fTrackID[t];
508 if(label < 0) break;
509 UInt_t c;
510 for(c=0; c<(MaxTrack-1); c++)
917e711b 511 if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
0bd0c1ef 512 break;
513 // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
917e711b 514 etaclust[etaindex].fMcLabels[c] = label;
0bd0c1ef 515 }
516 }
517#endif
518 continue;
519 }
520 else
521 {
917e711b 522 if(pad <= (etaclust[etaindex].fEndPad + 1))
0bd0c1ef 523 {
917e711b 524 etaclust[etaindex].fEndPad = pad;
0bd0c1ef 525#ifdef do_mc
526 if(fDoMC)
527 {
528 for(Int_t t=0; t<3; t++)
529 {
530 Int_t label = digPt[j].fTrackID[t];
531 if(label < 0) break;
532 UInt_t c;
533 for(c=0; c<(MaxTrack-1); c++)
917e711b 534 if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
0bd0c1ef 535 break;
536 // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
917e711b 537 etaclust[etaindex].fMcLabels[c] = label;
0bd0c1ef 538 }
539 }
540#endif
541 }
542 else
543 {
917e711b 544 Bool_t fillcluster = kTRUE;
545 if(fillcluster) {
0bd0c1ef 546
917e711b 547 UChar_t *nrows = fgRowCount[etaindex];
548 UChar_t *ngaps = fgGapCount[etaindex];
549 UChar_t *currentrow = fgCurrentRowCount[etaindex];
550 UChar_t *lastrow = fgTrackLastRow;
0bd0c1ef 551
552 //Do the transformation:
917e711b 553 Float_t starty = etaclust[etaindex].fStartY;
554 if(etaclust[etaindex].fStartPad == 0)
555 starty -= 0.0*padpitch;
556 Float_t r1 = x2 + starty*starty;
557 Float_t xoverr1 = x/r1;
558 Float_t startyoverr1 = starty/r1;
559 Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
560 if(etaclust[etaindex].fEndPad == npads)
561 endy += 0.0*padpitch;
562 Float_t r2 = x2 + endy*endy;
563 Float_t xoverr2 = x/r2;
564 Float_t endyoverr2 = endy/r2;
565 Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
566 Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
567 Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
568 Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
569
e81040b7 570 Float_t alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
917e711b 571 Float_t deltaalpha1 = b1*histbin/xbin;
572 if(b1<0)
e81040b7 573 alpha1 += deltaalpha1;
574 Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
917e711b 575 Float_t deltaalpha2 = b2*histbin/xbin;
576 if(b2>=0)
e81040b7 577 alpha2 += deltaalpha2;
0bd0c1ef 578
579 //Fill the histogram along the phirange
e81040b7 580 for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
0bd0c1ef 581 {
e81040b7 582 Int_t binx1 = 1 + (Int_t)alpha1;
917e711b 583 if(binx1>lastbinx) continue;
584 if(binx1<firstbinx) binx1 = firstbinx;
e81040b7 585 Int_t binx2 = 1 + (Int_t)alpha2;
917e711b 586 if(binx2<firstbinx) continue;
587 if(binx2>lastbinx) binx2 = lastbinx;
de3c3890 588#ifdef do_mc
589 if(binx2<binx1) {
590 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
917e711b 591 <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
de3c3890 592 }
593#endif
917e711b 594 Int_t tempbin = b*nbinx;
595 UChar_t *nrows2 = nrows + tempbin;
596 UChar_t *ngaps2 = ngaps + tempbin;
597 UChar_t *currentrow2 = currentrow + tempbin;
598 UChar_t *lastrow2 = lastrow + tempbin;
0bd0c1ef 599 for(Int_t bin=binx1;bin<=binx2;bin++)
600 {
de3c3890 601 if(ngaps2[bin] < MAX_N_GAPS) {
602 if(i > (currentrow2[bin] + MAX_GAP_SIZE) && i < lastrow2[bin]) {
603 ngaps2[bin] = MAX_N_GAPS;
604 continue;
605 }
606 if(i > currentrow2[bin] && i < lastrow2[bin])
0bd0c1ef 607 {
de3c3890 608 nrows2[bin]++;
609 if(i > (currentrow2[bin] + 1))
610 ngaps2[bin]++;
611 currentrow2[bin]=i;
0bd0c1ef 612 }
613#ifdef do_mc
614 if(fDoMC)
615 {
616 for(UInt_t t=0;t<(MaxTrack-1); t++)
617 {
917e711b 618 Int_t label = etaclust[etaindex].fMcLabels[t];
0bd0c1ef 619 if(label == 0) break;
620 UInt_t c;
917e711b 621 Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
0bd0c1ef 622 for(c=0; c<(MaxTrack-1); c++)
917e711b 623 if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
0bd0c1ef 624 break;
625 // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
917e711b 626 fgTrackID[etaindex][tempbin2].fLabel[c] = label;
627 if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
628 fgTrackID[etaindex][tempbin2].fNHits[c]++;
629 fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
0bd0c1ef 630 }
631 }
632 }
633#endif
634 }
0bd0c1ef 635 }
636 }
637 //End of the transformation
638
0bd0c1ef 639 }
0bd0c1ef 640
917e711b 641 etaclust[etaindex].fStartPad = pad;
642 etaclust[etaindex].fEndPad = pad;
643 etaclust[etaindex].fStartY = y - padpitch/2.0;
0bd0c1ef 644
645#ifdef do_mc
646 if(fDoMC)
647 {
917e711b 648 memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
0bd0c1ef 649 for(Int_t t=0; t<3; t++)
650 {
651 Int_t label = digPt[j].fTrackID[t];
652 if(label < 0) break;
653 UInt_t c;
654 for(c=0; c<(MaxTrack-1); c++)
917e711b 655 if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
0bd0c1ef 656 break;
657 // if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
917e711b 658 etaclust[etaindex].fMcLabels[c] = label;
0bd0c1ef 659 }
660 }
661#endif
662 }
663 }
917e711b 664 lastpad = pad;
665 lastetaindex = etaindex;
0bd0c1ef 666 }
667 //Write remaining clusters
917e711b 668 for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
0bd0c1ef 669 {
670 //Check for empty row
917e711b 671 if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
0bd0c1ef 672
917e711b 673 UChar_t *nrows = fgRowCount[etaindex];
674 UChar_t *ngaps = fgGapCount[etaindex];
675 UChar_t *currentrow = fgCurrentRowCount[etaindex];
676 UChar_t *lastrow = fgTrackLastRow;
0bd0c1ef 677
678 //Do the transformation:
917e711b 679 Float_t starty = etaclust[etaindex].fStartY;
680 if(etaclust[etaindex].fStartPad == 0)
681 starty -= 0.0*padpitch;
682 Float_t r1 = x2 + starty*starty;
683 Float_t xoverr1 = x/r1;
684 Float_t startyoverr1 = starty/r1;
685 Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
686 if(etaclust[etaindex].fEndPad == npads)
687 endy += 0.0*padpitch;
688 Float_t r2 = x2 + endy*endy;
689 Float_t xoverr2 = x/r2;
690 Float_t endyoverr2 = endy/r2;
691 Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
692 Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
693 Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
694 Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
695
e81040b7 696 Float_t alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
917e711b 697 Float_t deltaalpha1 = b1*histbin/xbin;
698 if(b1<0)
e81040b7 699 alpha1 += deltaalpha1;
700 Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
917e711b 701 Float_t deltaalpha2 = b2*histbin/xbin;
702 if(b2>=0)
e81040b7 703 alpha2 += deltaalpha2;
0bd0c1ef 704
705 //Fill the histogram along the phirange
e81040b7 706 for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
0bd0c1ef 707 {
e81040b7 708 Int_t binx1 = 1 + (Int_t)alpha1;
917e711b 709 if(binx1>lastbinx) continue;
710 if(binx1<firstbinx) binx1 = firstbinx;
e81040b7 711 Int_t binx2 = 1 + (Int_t)alpha2;
917e711b 712 if(binx2<firstbinx) continue;
713 if(binx2>lastbinx) binx2 = lastbinx;
de3c3890 714#ifdef do_mc
715 if(binx2<binx1) {
716 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
917e711b 717 <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
de3c3890 718 }
719#endif
917e711b 720 Int_t tempbin = b*nbinx;
721 UChar_t *nrows2 = nrows + tempbin;
722 UChar_t *ngaps2 = ngaps + tempbin;
723 UChar_t *currentrow2 = currentrow + tempbin;
724 UChar_t *lastrow2 = lastrow + tempbin;
0bd0c1ef 725 for(Int_t bin=binx1;bin<=binx2;bin++)
726 {
de3c3890 727 if(ngaps2[bin] < MAX_N_GAPS) {
728 if(i > (currentrow2[bin] + MAX_GAP_SIZE) && i < lastrow2[bin]) {
729 ngaps2[bin] = MAX_N_GAPS;
730 continue;
731 }
732 if(i > currentrow2[bin] && i < lastrow2[bin])
0bd0c1ef 733 {
de3c3890 734 nrows2[bin]++;
735 if(i > (currentrow2[bin] + 1))
736 ngaps2[bin]++;
737 currentrow2[bin]=i;
0bd0c1ef 738 }
739#ifdef do_mc
740 if(fDoMC)
741 {
742 for(UInt_t t=0;t<(MaxTrack-1); t++)
743 {
917e711b 744 Int_t label = etaclust[etaindex].fMcLabels[t];
0bd0c1ef 745 if(label == 0) break;
746 UInt_t c;
917e711b 747 Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
0bd0c1ef 748 for(c=0; c<(MaxTrack-1); c++)
917e711b 749 if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
0bd0c1ef 750 break;
751 // if(c == MaxTrack-1) cout<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
917e711b 752 fgTrackID[etaindex][tempbin2].fLabel[c] = label;
753 if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
754 fgTrackID[etaindex][tempbin2].fNHits[c]++;
755 fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
0bd0c1ef 756 }
757 }
758 }
759#endif
760 }
0bd0c1ef 761 }
762 }
763 //End of the transformation
764
0bd0c1ef 765 }
766
767 //Move the data pointer to the next padrow:
768 AliL3MemHandler::UpdateRowPointer(tempPt);
769 }
770
917e711b 771 delete [] etaclust;
0bd0c1ef 772}
773
e81040b7 774Int_t AliL3HoughTransformerRow::GetTrackID(Int_t etaindex,Double_t alpha1,Double_t alpha2) const
0bd0c1ef 775{
917e711b 776 // Returns the MC label for a given peak found in the Hough space
0bd0c1ef 777 if(!fDoMC)
778 {
de3c3890 779 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
780 <<"Flag switched off"<<ENDLOG;
0bd0c1ef 781 return -1;
782 }
783
784#ifdef do_mc
917e711b 785 if(etaindex < 0 || etaindex > GetNEtaSegments())
0bd0c1ef 786 {
de3c3890 787 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
917e711b 788 <<"Wrong etaindex "<<etaindex<<ENDLOG;
0bd0c1ef 789 return -1;
790 }
917e711b 791 AliL3Histogram *hist = fParamSpace[etaindex];
e81040b7 792 Int_t bin = hist->FindLabelBin(alpha1,alpha2);
e81ffe36 793 if(bin==-1) {
794 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
e81040b7 795 <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
e81ffe36 796 return -1;
797 }
0bd0c1ef 798 Int_t label=-1;
799 Int_t max=0;
800 for(UInt_t i=0; i<(MaxTrack-1); i++)
801 {
917e711b 802 Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
0bd0c1ef 803 if(nhits == 0) break;
804 if(nhits > max)
805 {
806 max = nhits;
917e711b 807 label = fgTrackID[etaindex][bin].fLabel[i];
0bd0c1ef 808 }
809 }
810 Int_t label2=-1;
811 Int_t max2=0;
812 for(UInt_t i=0; i<(MaxTrack-1); i++)
813 {
917e711b 814 Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
0bd0c1ef 815 if(nhits == 0) break;
816 if(nhits > max2)
817 {
917e711b 818 if(fgTrackID[etaindex][bin].fLabel[i]!=label) {
0bd0c1ef 819 max2 = nhits;
917e711b 820 label2 = fgTrackID[etaindex][bin].fLabel[i];
0bd0c1ef 821 }
822 }
823 }
917e711b 824 if(max2 !=0 ) {
825 LOG(AliL3Log::kDebug,"AliL3HoughTransformerRow::GetTrackID()","")
826 <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fgTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fgTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
827 }
0bd0c1ef 828 return label;
829#endif
de3c3890 830 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
831 <<"Compile with do_mc flag!"<<ENDLOG;
0bd0c1ef 832 return -1;
833}
834
917e711b 835UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t etaindex)
0bd0c1ef 836{
917e711b 837 return fgRowCount[etaindex];
0bd0c1ef 838}
839
917e711b 840UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t etaindex)
0bd0c1ef 841{
917e711b 842 return fgGapCount[etaindex];
0bd0c1ef 843}
de3c3890 844
917e711b 845UChar_t *AliL3HoughTransformerRow::GetCurrentRowCount(Int_t etaindex)
de3c3890 846{
917e711b 847 return fgCurrentRowCount[etaindex];
de3c3890 848}
849
850UChar_t *AliL3HoughTransformerRow::GetTrackNRows()
851{
917e711b 852 return fgTrackNRows;
de3c3890 853}
854
855
856UChar_t *AliL3HoughTransformerRow::GetTrackFirstRow()
857{
917e711b 858 return fgTrackFirstRow;
de3c3890 859}
860
861UChar_t *AliL3HoughTransformerRow::GetTrackLastRow()
862{
917e711b 863 return fgTrackLastRow;
de3c3890 864}