2 // origin hough/AliL3HoughTransformerRow.cxx,v 1.21 Tue Mar 28 18:05:12 2006 UTC by alibrary
4 /** @file AliHLTTPCHoughTransformerRow.cxx
5 @author Cvetan Cheshkov <mailto:cvetan.cheshkov@cern.ch>
7 @brief Implementation of fast HLT TPC hough transform tracking. */
9 //_____________________________________________________________
10 // AliHLTTPCHoughTransformerRow
12 // Impelementation of the so called "TPC rows Hough transformation" class
14 // Transforms the TPC data into the hough space and counts the missed TPC
15 // rows corresponding to each track cnadidate - hough space bin
17 #include "AliHLTStdIncludes.h"
19 #include "AliHLTTPCLogging.h"
20 #include "AliHLTTPCMemHandler.h"
21 #include "AliHLTTPCTransform.h"
22 #include "AliHLTTPCDigitData.h"
23 #include "AliHLTTPCHistogramAdaptive.h"
24 #include "AliHLTTPCHoughTrack.h"
25 #include "AliHLTTPCHoughTransformerRow.h"
26 #include "AliTPCRawStream.h"
32 /** ROOT macro for the implementation of ROOT specific class methods */
33 ClassImp(AliHLTTPCHoughTransformerRow);
35 Float_t AliHLTTPCHoughTransformerRow::fgBeta1 = 1.0/AliHLTTPCTransform::Row2X(84);
36 Float_t AliHLTTPCHoughTransformerRow::fgBeta2 = 1.0/(AliHLTTPCTransform::Row2X(158)*(1.0+tan(AliHLTTPCTransform::Pi()*10/180)*tan(AliHLTTPCTransform::Pi()*10/180)));
37 Float_t AliHLTTPCHoughTransformerRow::fgDAlpha = 0.22;
38 Float_t AliHLTTPCHoughTransformerRow::fgDEta = 0.40;
39 Double_t AliHLTTPCHoughTransformerRow::fgEtaCalcParam1 = 1.0289;
40 Double_t AliHLTTPCHoughTransformerRow::fgEtaCalcParam2 = 0.15192;
41 Double_t AliHLTTPCHoughTransformerRow::fgEtaCalcParam3 = 1./(32.*600.*600.);
43 AliHLTTPCHoughTransformerRow::AliHLTTPCHoughTransformerRow()
70 AliHLTTPCHoughTransformerRow::AliHLTTPCHoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliHLTTPCHoughTransformer(slice,patch,netasegments,zvertex)
98 AliHLTTPCHoughTransformerRow::~AliHLTTPCHoughTransformerRow()
101 if(fLastTransformer) return;
106 for(Int_t i=0; i<GetNEtaSegments(); i++)
108 if(!fTrackID[i]) continue;
118 for(Int_t i=0; i<GetNEtaSegments(); i++)
120 if(!fGapCount[i]) continue;
121 delete [] fGapCount[i];
128 for(Int_t i=0; i<GetNEtaSegments(); i++)
130 if(fCurrentRowCount[i])
131 delete [] fCurrentRowCount[i];
133 delete [] fCurrentRowCount;
134 fCurrentRowCount = 0;
138 delete [] fTrackNRows;
143 delete [] fTrackFirstRow;
148 delete [] fTrackLastRow;
153 delete [] fInitialGapCount;
154 fInitialGapCount = 0;
158 for(Int_t i=0; i<GetNEtaSegments(); i++)
160 if(!fPrevBin[i]) continue;
161 delete [] fPrevBin[i];
168 for(Int_t i=0; i<GetNEtaSegments(); i++)
170 if(!fNextBin[i]) continue;
171 delete [] fNextBin[i];
178 for(Int_t i=0; i<GetNEtaSegments(); i++)
180 if(!fNextRow[i]) continue;
181 delete [] fNextRow[i];
188 for(Int_t i = AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
190 if(!fStartPadParams[i]) continue;
191 delete [] fStartPadParams[i];
193 delete [] fStartPadParams;
198 for(Int_t i = AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
200 if(!fEndPadParams[i]) continue;
201 delete [] fEndPadParams[i];
203 delete [] fEndPadParams;
208 for(Int_t i = AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
210 if(!fLUTr[i]) continue;
218 delete[] fLUTforwardZ;
223 delete[] fLUTbackwardZ;
228 void AliHLTTPCHoughTransformerRow::DeleteHistograms()
233 for(Int_t i=0; i<GetNEtaSegments(); i++)
235 if(!fParamSpace[i]) continue;
236 delete fParamSpace[i];
238 delete [] fParamSpace;
241 void AliHLTTPCHoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
242 Int_t nybin,Float_t ymin,Float_t ymax)
244 //Create the histograms (parameter space)
247 //xmin xmax ymin ymax = histogram limits in X and Y
248 if(fLastTransformer) {
249 SetTransformerArrays((AliHLTTPCHoughTransformerRow *)fLastTransformer);
252 fParamSpace = new AliHLTTPCHistogram*[GetNEtaSegments()];
254 Char_t histname[256];
255 for(Int_t i=0; i<GetNEtaSegments(); i++)
257 sprintf(histname,"paramspace_%d",i);
258 fParamSpace[i] = new AliHLTTPCHistogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
262 AliHLTTPCHistogram *hist = fParamSpace[0];
263 Int_t ncellsx = (hist->GetNbinsX()+3)/2;
264 Int_t ncellsy = (hist->GetNbinsY()+3)/2;
265 Int_t ncells = ncellsx*ncellsy;
268 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
269 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliHLTTrackIndex)<<" bytes to fTrackID"<<ENDLOG;
270 fTrackID = new AliHLTTrackIndex*[GetNEtaSegments()];
271 for(Int_t i=0; i<GetNEtaSegments(); i++)
272 fTrackID[i] = new AliHLTTrackIndex[ncells];
276 AliHLTTPCHistogram *hist = fParamSpace[0];
277 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
280 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
281 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<ENDLOG;
282 fGapCount = new UChar_t*[GetNEtaSegments()];
283 for(Int_t i=0; i<GetNEtaSegments(); i++)
284 fGapCount[i] = new UChar_t[ncells];
286 if(!fCurrentRowCount)
288 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
289 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<ENDLOG;
290 fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
291 for(Int_t i=0; i<GetNEtaSegments(); i++)
292 fCurrentRowCount[i] = new UChar_t[ncells];
296 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
297 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fPrevBin"<<ENDLOG;
298 fPrevBin = new UChar_t*[GetNEtaSegments()];
299 for(Int_t i=0; i<GetNEtaSegments(); i++)
300 fPrevBin[i] = new UChar_t[ncells];
304 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
305 <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fNextBin"<<ENDLOG;
306 fNextBin = new UChar_t*[GetNEtaSegments()];
307 for(Int_t i=0; i<GetNEtaSegments(); i++)
308 fNextBin[i] = new UChar_t[ncells];
310 Int_t ncellsy = hist->GetNbinsY()+2;
313 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
314 <<"Transformer: Allocating "<<GetNEtaSegments()*ncellsy*sizeof(UChar_t)<<" bytes to fNextRow"<<ENDLOG;
315 fNextRow = new UChar_t*[GetNEtaSegments()];
316 for(Int_t i=0; i<GetNEtaSegments(); i++)
317 fNextRow[i] = new UChar_t[ncellsy];
322 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
323 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackNRows"<<ENDLOG;
324 fTrackNRows = new UChar_t[ncells];
325 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
326 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackFirstRow"<<ENDLOG;
327 fTrackFirstRow = new UChar_t[ncells];
328 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
329 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackLastRow"<<ENDLOG;
330 fTrackLastRow = new UChar_t[ncells];
331 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
332 <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fInitialGapCount"<<ENDLOG;
333 fInitialGapCount = new UChar_t[ncells];
336 AliHLTTPCHoughTrack track;
337 Int_t xmin = hist->GetFirstXbin();
338 Int_t xmax = hist->GetLastXbin();
339 Int_t xmiddle = (hist->GetNbinsX()+1)/2;
340 Int_t ymin = hist->GetFirstYbin();
341 Int_t ymax = hist->GetLastYbin();
342 Int_t nxbins = hist->GetNbinsX()+2;
343 Int_t nxgrid = (hist->GetNbinsX()+3)/2+1;
344 Int_t nygrid = hist->GetNbinsY()+3;
346 AliHLTTrackLength *tracklength = new AliHLTTrackLength[nxgrid*nygrid];
347 memset(tracklength,0,nxgrid*nygrid*sizeof(AliHLTTrackLength));
349 for(Int_t ybin=ymin-1; ybin<=(ymax+1); ybin++)
351 for(Int_t xbin=xmin-1; xbin<=xmiddle; xbin++)
353 fTrackNRows[xbin + ybin*nxbins] = 255;
354 for(Int_t deltay = 0; deltay <= 1; deltay++) {
355 for(Int_t deltax = 0; deltax <= 1; deltax++) {
357 AliHLTTrackLength *curtracklength = &tracklength[(xbin + deltax) + (ybin + deltay)*nxgrid];
358 UInt_t maxfirstrow = 0;
359 UInt_t maxlastrow = 0;
360 Float_t maxtrackpt = 0;
361 if(curtracklength->fIsFilled) {
362 maxfirstrow = curtracklength->fFirstRow;
363 maxlastrow = curtracklength->fLastRow;
364 maxtrackpt = curtracklength->fTrackPt;
367 Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)(2*deltax-1));
368 Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)(2*deltay-1));
370 Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
371 Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
372 track.SetTrackParameters(kappa,psi,1);
373 maxtrackpt = track.GetPt();
374 if(maxtrackpt < 0.9*0.1*AliHLTTPCTransform::GetSolenoidField())
376 maxfirstrow = maxlastrow = 0;
377 curtracklength->fIsFilled = kTRUE;
378 curtracklength->fFirstRow = maxfirstrow;
379 curtracklength->fLastRow = maxlastrow;
380 curtracklength->fTrackPt = maxtrackpt;
384 Bool_t firstrow = kFALSE;
385 UInt_t curfirstrow = 0;
386 UInt_t curlastrow = 0;
388 Double_t centerx = track.GetCenterX();
389 Double_t centery = track.GetCenterY();
390 Double_t radius = track.GetRadius();
392 for(Int_t j=AliHLTTPCTransform::GetFirstRow(0); j<=AliHLTTPCTransform::GetLastRow(5); j++)
395 // if(!track.GetCrossingPoint(j,hit)) continue;
396 hit[0] = AliHLTTPCTransform::Row2X(j);
397 Double_t aa = (hit[0] - centerx)*(hit[0] - centerx);
398 Double_t r2 = radius*radius;
402 Double_t aa2 = sqrt(r2 - aa);
403 Double_t y1 = centery + aa2;
404 Double_t y2 = centery - aa2;
406 if(fabs(y2) < fabs(y1)) hit[1] = y2;
410 AliHLTTPCTransform::LocHLT2Raw(hit,0,j);
412 if(hit[1]>=0 && hit[1]<AliHLTTPCTransform::GetNPads(j))
423 if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
424 maxfirstrow = curfirstrow;
425 maxlastrow = curlastrow;
430 if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
431 maxfirstrow = curfirstrow;
432 maxlastrow = curlastrow;
435 curtracklength->fIsFilled = kTRUE;
436 curtracklength->fFirstRow = maxfirstrow;
437 curtracklength->fLastRow = maxlastrow;
438 curtracklength->fTrackPt = maxtrackpt;
441 if((maxlastrow-maxfirstrow) < fTrackNRows[xbin + ybin*nxbins]) {
442 fTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
443 fInitialGapCount[xbin + ybin*nxbins] = 1;
444 if((maxlastrow-maxfirstrow+1)<=MIN_TRACK_LENGTH)
445 fInitialGapCount[xbin + ybin*nxbins] = MAX_N_GAPS+1;
446 if(maxtrackpt < 0.9*0.1*AliHLTTPCTransform::GetSolenoidField())
447 fInitialGapCount[xbin + ybin*nxbins] = MAX_N_GAPS;
448 fTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
449 fTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
451 Int_t xmirror = xmax - xbin + 1;
452 Int_t ymirror = ymax - ybin + 1;
453 fTrackNRows[xmirror + ymirror*nxbins] = fTrackNRows[xbin + ybin*nxbins];
454 fInitialGapCount[xmirror + ymirror*nxbins] = fInitialGapCount[xbin + ybin*nxbins];
455 fTrackFirstRow[xmirror + ymirror*nxbins] = fTrackFirstRow[xbin + ybin*nxbins];
456 fTrackLastRow[xmirror + ymirror*nxbins] = fTrackLastRow[xbin + ybin*nxbins];
460 // cout<<" fTrackNRows "<<(Int_t)fInitialGapCount[xbin + ybin*nxbins]<<" "<<xbin<<" "<<ybin<<" "<<(Int_t)fTrackNRows[xbin + ybin*nxbins]<<" "<<(Int_t)fTrackFirstRow[xbin + ybin*nxbins]<<" "<<(Int_t)fTrackLastRow[xbin + ybin*nxbins]<<" "<<endl;
463 delete [] tracklength;
468 Int_t nrows = AliHLTTPCTransform::GetLastRow(5) - AliHLTTPCTransform::GetFirstRow(0) + 1;
469 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
470 <<"Transformer: Allocating about "<<nrows*100*sizeof(AliHLTPadHoughParams)<<" bytes to fStartPadParams"<<ENDLOG;
471 fStartPadParams = new AliHLTPadHoughParams*[nrows];
472 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
473 <<"Transformer: Allocating about "<<nrows*100*sizeof(AliHLTPadHoughParams)<<" bytes to fEndPadParams"<<ENDLOG;
474 fEndPadParams = new AliHLTPadHoughParams*[nrows];
475 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
476 <<"Transformer: Allocating about "<<nrows*100*sizeof(Float_t)<<" bytes to fLUTr"<<ENDLOG;
477 fLUTr = new Float_t*[nrows];
479 Float_t beta1 = fgBeta1;
480 Float_t beta2 = fgBeta2;
481 Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
482 Float_t ymin = hist->GetYmin();
483 Float_t histbin = hist->GetBinWidthY();
484 Float_t xmin = hist->GetXmin();
485 Float_t xmax = hist->GetXmax();
486 Float_t xbin = (xmax-xmin)/hist->GetNbinsX();
487 Int_t firstbinx = hist->GetFirstXbin();
488 Int_t lastbinx = hist->GetLastXbin();
489 Int_t nbinx = hist->GetNbinsX()+2;
490 Int_t firstbin = hist->GetFirstYbin();
491 Int_t lastbin = hist->GetLastYbin();
492 for(Int_t i=AliHLTTPCTransform::GetFirstRow(0); i<=AliHLTTPCTransform::GetLastRow(5); i++)
494 Int_t npads = AliHLTTPCTransform::GetNPads(i);
495 Int_t ipatch = AliHLTTPCTransform::GetPatch(i);
496 Double_t padpitch = AliHLTTPCTransform::GetPadPitchWidth(ipatch);
497 Float_t x = AliHLTTPCTransform::Row2X(i);
500 fStartPadParams[i] = new AliHLTPadHoughParams[npads];
501 fEndPadParams[i] = new AliHLTPadHoughParams[npads];
502 fLUTr[i] = new Float_t[npads];
503 for(Int_t pad=0; pad<npads; pad++)
505 Float_t y = (pad-0.5*(npads-1))*padpitch;
506 fLUTr[i][pad] = sqrt(x2 + y*y);
507 Float_t starty = (pad-0.5*npads)*padpitch;
508 Float_t r1 = x2 + starty*starty;
509 Float_t xoverr1 = x/r1;
510 Float_t startyoverr1 = starty/r1;
511 Float_t endy = (pad-0.5*(npads-2))*padpitch;
512 Float_t r2 = x2 + endy*endy;
513 Float_t xoverr2 = x/r2;
514 Float_t endyoverr2 = endy/r2;
515 Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
516 Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
517 Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
518 Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
520 Float_t alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
521 Float_t deltaalpha1 = b1*histbin/xbin;
523 alpha1 += deltaalpha1;
524 Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
525 Float_t deltaalpha2 = b2*histbin/xbin;
527 alpha2 += deltaalpha2;
529 fStartPadParams[i][pad].fAlpha = alpha1;
530 fStartPadParams[i][pad].fDeltaAlpha = deltaalpha1;
531 fEndPadParams[i][pad].fAlpha = alpha2;
532 fEndPadParams[i][pad].fDeltaAlpha = deltaalpha2;
534 //Find the first and last bin rows to be filled
535 Bool_t binfound1 = kFALSE;
536 Bool_t binfound2 = kFALSE;
537 Int_t firstbin1 = lastbin;
538 Int_t firstbin2 = lastbin;
539 Int_t lastbin1 = firstbin;
540 Int_t lastbin2 = firstbin;
541 for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
543 Int_t binx1 = 1 + (Int_t)alpha1;
544 if(binx1<=lastbinx) {
545 UChar_t initialgapcount;
547 initialgapcount = fInitialGapCount[binx1 + b*nbinx];
549 initialgapcount = fInitialGapCount[firstbinx + b*nbinx];
550 if(initialgapcount != MAX_N_GAPS) {
558 Int_t binx2 = 1 + (Int_t)alpha2;
559 if(binx2>=firstbin) {
560 UChar_t initialgapcount;
562 initialgapcount = fInitialGapCount[binx2 + b*nbinx];
564 initialgapcount = fInitialGapCount[lastbinx + b*nbinx];
565 if(initialgapcount != MAX_N_GAPS) {
574 fStartPadParams[i][pad].fFirstBin=firstbin1;
575 fStartPadParams[i][pad].fLastBin=lastbin1;
576 fEndPadParams[i][pad].fFirstBin=firstbin2;
577 fEndPadParams[i][pad].fLastBin=lastbin2;
582 //create lookup table for z of the digits
585 Int_t ntimebins = AliHLTTPCTransform::GetNTimeBins();
586 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
587 <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTforwardZ"<<ENDLOG;
588 fLUTforwardZ = new Float_t[ntimebins];
589 LOG(AliHLTTPCLog::kInformational,"AliHLTTPCHoughTransformerRow::CreateHistograms()","")
590 <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTbackwardZ"<<ENDLOG;
591 fLUTbackwardZ = new Float_t[ntimebins];
592 for(Int_t i=0; i<ntimebins; i++){
594 z=AliHLTTPCTransform::GetZFast(0,i,GetZVertex());
596 z=AliHLTTPCTransform::GetZFast(18,i,GetZVertex());
602 void AliHLTTPCHoughTransformerRow::Reset()
604 //Reset all the histograms. Should be done when processing new slice
605 if(fLastTransformer) return;
608 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformer::Reset","Histograms")
609 <<"No histograms to reset"<<ENDLOG;
613 for(Int_t i=0; i<GetNEtaSegments(); i++)
614 fParamSpace[i]->Reset();
618 AliHLTTPCHistogram *hist = fParamSpace[0];
619 Int_t ncellsx = (hist->GetNbinsX()+3)/2;
620 Int_t ncellsy = (hist->GetNbinsY()+3)/2;
621 Int_t ncells = ncellsx*ncellsy;
622 for(Int_t i=0; i<GetNEtaSegments(); i++)
623 memset(fTrackID[i],0,ncells*sizeof(AliHLTTrackIndex));
626 AliHLTTPCHistogram *hist = fParamSpace[0];
627 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
628 for(Int_t i=0; i<GetNEtaSegments(); i++)
630 memcpy(fGapCount[i],fInitialGapCount,ncells*sizeof(UChar_t));
631 memcpy(fCurrentRowCount[i],fTrackFirstRow,ncells*sizeof(UChar_t));
635 Int_t AliHLTTPCHoughTransformerRow::GetEtaIndex(Double_t eta) const
637 //Return the histogram index of the corresponding eta.
639 Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
640 Double_t index = (eta-GetEtaMin())/etaslice;
644 inline AliHLTTPCHistogram *AliHLTTPCHoughTransformerRow::GetHistogram(Int_t etaindex)
646 // Return a pointer to the histogram which contains etaindex eta slice
647 if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
649 if(!fParamSpace[etaindex])
651 return fParamSpace[etaindex];
654 Double_t AliHLTTPCHoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/) const
656 // Return eta calculated in the middle of the eta slice
657 Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
659 eta=(Double_t)((etaindex+0.5)*etaslice);
663 void AliHLTTPCHoughTransformerRow::TransformCircle()
665 // This method contains the hough transformation
666 // Depending on the option selected, it reads as an input
667 // either the preloaded array with digits or directly raw
668 // data stream (option fast_raw)
670 TransformCircleFromDigitArray();
671 else if(fTPCRawStream)
672 TransformCircleFromRawStream();
674 void AliHLTTPCHoughTransformerRow::TransformCircleFromDigitArray()
676 //Do the Hough Transform
678 //Load the parameters used by the fast calculation of eta
679 Double_t etaparam1 = GetEtaCalcParam1();
680 Double_t etaparam2 = GetEtaCalcParam2();
681 Double_t etaparam3 = GetEtaCalcParam3();
683 Int_t netasegments = GetNEtaSegments();
684 Double_t etamin = GetEtaMin();
685 Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
687 Int_t lowerthreshold = GetLowerThreshold();
689 //Assumes that all the etaslice histos are the same!
690 AliHLTTPCHistogram *h = fParamSpace[0];
691 Int_t firstbiny = h->GetFirstYbin();
692 Int_t firstbinx = h->GetFirstXbin();
693 Int_t lastbinx = h->GetLastXbin();
694 Int_t nbinx = h->GetNbinsX()+2;
697 Int_t lastetaindex=-1;
698 AliHLTEtaRow *etaclust = new AliHLTEtaRow[netasegments];
700 AliHLTTPCDigitRowData *tempPt = GetDataPointer();
703 LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTransformer::TransformCircle","Data")
704 <<"No input data "<<ENDLOG;
708 Int_t ipatch = GetPatch();
709 Int_t ilastpatch = GetLastPatch();
710 Int_t islice = GetSlice();
715 lutz = fLUTbackwardZ;
717 //Loop over the padrows:
718 for(UChar_t i=AliHLTTPCTransform::GetFirstRow(ipatch); i<=AliHLTTPCTransform::GetLastRow(ipatch); i++)
721 //Flush eta clusters array
722 memset(etaclust,0,netasegments*sizeof(AliHLTEtaRow));
726 //Get the data on this padrow:
727 AliHLTTPCDigitData *digPt = tempPt->fDigitData;
728 if((Int_t)i != (Int_t)tempPt->fRow)
730 LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTransformerRow::TransformCircle","Data")
731 <<"AliHLTTPCHoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<ENDLOG;
734 // cout<<" Starting row "<<i<<endl;
735 //Loop over the data on this padrow:
736 for(UInt_t j=0; j<tempPt->fNDigit; j++)
738 UShort_t charge = digPt[j].fCharge;
739 if((Int_t)charge <= lowerthreshold)
741 UChar_t pad = digPt[j].fPad;
742 UShort_t time = digPt[j].fTime;
746 radius = fLUTr[(Int_t)i][(Int_t)pad];
750 Float_t z = lutz[(Int_t)time];
751 Double_t radiuscorr = radius*(1.+etaparam3*radius*radius);
752 Double_t zovr = z/radiuscorr;
753 Double_t eta = (etaparam1-etaparam2*fabs(zovr))*zovr;
754 //Get the corresponding index, which determines which histogram to fill:
755 Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
758 if(etaindex == lastetaindex) continue;
760 // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
762 if(etaindex < 0 || etaindex >= netasegments)
765 if(!etaclust[etaindex].fIsFound)
767 etaclust[etaindex].fStartPad = pad;
768 etaclust[etaindex].fEndPad = pad;
769 etaclust[etaindex].fIsFound = 1;
771 FillClusterMCLabels(digPt[j],&etaclust[etaindex]);
777 if(pad <= (etaclust[etaindex].fEndPad + 1))
779 etaclust[etaindex].fEndPad = pad;
781 FillClusterMCLabels(digPt[j],&etaclust[etaindex]);
786 FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
788 etaclust[etaindex].fStartPad = pad;
789 etaclust[etaindex].fEndPad = pad;
792 memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
793 FillClusterMCLabels(digPt[j],&etaclust[etaindex]);
798 lastetaindex = etaindex;
800 //Write remaining clusters
801 for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
803 //Check for empty row
804 if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
806 FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
809 //Move the data pointer to the next padrow:
810 AliHLTTPCMemHandler::UpdateRowPointer(tempPt);
816 void AliHLTTPCHoughTransformerRow::TransformCircleFromRawStream()
818 //Do the Hough Transform
820 //Load the parameters used by the fast calculation of eta
821 Double_t etaparam1 = GetEtaCalcParam1();
822 Double_t etaparam2 = GetEtaCalcParam2();
823 Double_t etaparam3 = GetEtaCalcParam3();
825 Int_t netasegments = GetNEtaSegments();
826 Double_t etamin = GetEtaMin();
827 Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
829 Int_t lowerthreshold = GetLowerThreshold();
831 //Assumes that all the etaslice histos are the same!
832 AliHLTTPCHistogram *h = fParamSpace[0];
833 Int_t firstbiny = h->GetFirstYbin();
834 Int_t firstbinx = h->GetFirstXbin();
835 Int_t lastbinx = h->GetLastXbin();
836 Int_t nbinx = h->GetNbinsX()+2;
838 Int_t lastetaindex = -1;
839 AliHLTEtaRow *etaclust = new AliHLTEtaRow[netasegments];
843 LOG(AliHLTTPCLog::kError,"AliHLTTPCHoughTransformer::TransformCircle","Data")
844 <<"No input data "<<ENDLOG;
848 Int_t ipatch = GetPatch();
849 UChar_t rowmin = AliHLTTPCTransform::GetFirstRowOnDDL(ipatch);
850 UChar_t rowmax = AliHLTTPCTransform::GetLastRowOnDDL(ipatch);
851 // Int_t ntimebins = AliHLTTPCTransform::GetNTimeBins();
852 Int_t ilastpatch = GetLastPatch();
853 Int_t islice = GetSlice();
858 lutz = fLUTbackwardZ;
860 //Flush eta clusters array
861 memset(etaclust,0,netasegments*sizeof(AliHLTEtaRow));
868 //Loop over the rawdata:
869 while (fTPCRawStream->Next()) {
871 if(fTPCRawStream->IsNewSector() || fTPCRawStream->IsNewRow()) {
873 //Write remaining clusters
874 for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
876 //Check for empty row
877 if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
879 FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
882 Int_t sector=fTPCRawStream->GetSector();
883 Int_t row=fTPCRawStream->GetRow();
885 AliHLTTPCTransform::Sector2Slice(slice,srow,sector,row);
887 LOG(AliHLTTPCLog::kError,"AliHLTDDLDataFileHandler::DDLDigits2Memory","Slice")
888 <<AliHLTTPCLog::kDec<<"Found slice "<<slice<<", expected "<<islice<<ENDLOG;
893 npads = AliHLTTPCTransform::GetNPads(srow)-1;
895 //Flush eta clusters array
896 memset(etaclust,0,netasegments*sizeof(AliHLTEtaRow));
902 if((i<rowmin)||(i>rowmax))continue;
904 // cout<<" Starting row "<<(UInt_t)i<<endl;
905 //Loop over the data on this padrow:
906 if(fTPCRawStream->IsNewRow() || fTPCRawStream->IsNewPad()) {
907 pad=fTPCRawStream->GetPad();
909 if((pad<0)||(pad>=(npads+1))){
910 LOG(AliHLTTPCLog::kError,"AliHLTDDLDataFileHandler::DDLDigits2Memory","Pad")
911 <<AliHLTTPCLog::kDec<<"Pad value out of bounds "<<pad<<" "
916 radius = fLUTr[(Int_t)i][(Int_t)pad];
920 UShort_t time=fTPCRawStream->GetTime();
922 if((time<0)||(time>=ntimebins)){
923 LOG(AliHLTTPCLog::kError,"AliHLTDDLDataFileHandler::DDLDigits2Memory","Time")
924 <<AliHLTTPCLog::kDec<<"Time out of bounds "<<time<<" "
925 <<AliHLTTPCTransform::GetNTimeBins()<<ENDLOG;
930 if(fTPCRawStream->GetSignal() <= lowerthreshold)
933 Float_t z = lutz[(Int_t)time];
934 Double_t radiuscorr = radius*(1.+etaparam3*radius*radius);
935 Double_t zovr = z/radiuscorr;
936 Double_t eta = (etaparam1-etaparam2*fabs(zovr))*zovr;
937 //Get the corresponding index, which determines which histogram to fill:
938 Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
941 if(etaindex == lastetaindex) continue;
943 // cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
945 if(etaindex < 0 || etaindex >= netasegments)
948 if(!etaclust[etaindex].fIsFound)
950 etaclust[etaindex].fStartPad = pad;
951 etaclust[etaindex].fEndPad = pad;
952 etaclust[etaindex].fIsFound = 1;
957 if(pad <= (etaclust[etaindex].fEndPad + 1))
959 etaclust[etaindex].fEndPad = pad;
963 FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
965 etaclust[etaindex].fStartPad = pad;
966 etaclust[etaindex].fEndPad = pad;
970 lastetaindex = etaindex;
973 //Write remaining clusters
974 for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
976 //Check for empty row
977 if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
979 FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
986 Int_t AliHLTTPCHoughTransformerRow::GetTrackID(Int_t /*etaindex*/,Double_t /*alpha1*/,Double_t /*alpha2*/) const
988 // Does nothing if do_mc undefined
989 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackID","Data")
990 <<"Flag switched off"<<ENDLOG;
993 Int_t AliHLTTPCHoughTransformerRow::GetTrackID(Int_t etaindex,Double_t alpha1,Double_t alpha2) const
995 // Returns the MC label for a given peak found in the Hough space
996 if(etaindex < 0 || etaindex > GetNEtaSegments())
998 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackID","Data")
999 <<"Wrong etaindex "<<etaindex<<ENDLOG;
1002 AliHLTTPCHistogram *hist = fParamSpace[etaindex];
1003 Int_t bin = hist->FindLabelBin(alpha1,alpha2);
1005 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackID()","")
1006 <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
1011 for(UInt_t i=0; i<(MaxTrack-1); i++)
1013 Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
1014 if(nhits == 0) break;
1018 label = fTrackID[etaindex][bin].fLabel[i];
1023 for(UInt_t i=0; i<(MaxTrack-1); i++)
1025 Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
1026 if(nhits == 0) break;
1029 if(fTrackID[etaindex][bin].fLabel[i]!=label) {
1031 label2 = fTrackID[etaindex][bin].fLabel[i];
1036 LOG(AliHLTTPCLog::kDebug,"AliHLTTPCHoughTransformerRow::GetTrackID()","")
1037 <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
1043 Int_t AliHLTTPCHoughTransformerRow::GetTrackLength(Double_t alpha1,Double_t alpha2,Int_t *rows) const
1045 // Returns the track length for a given peak found in the Hough space
1047 AliHLTTPCHistogram *hist = fParamSpace[0];
1048 Int_t bin = hist->FindBin(alpha1,alpha2);
1050 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::GetTrackLength()","")
1051 <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
1054 rows[0] = fTrackFirstRow[bin];
1055 rows[1] = fTrackLastRow[bin];
1060 inline void AliHLTTPCHoughTransformerRow::FillClusterRow(UChar_t i,Int_t binx1,Int_t binx2,UChar_t *ngaps2,UChar_t *currentrow2,UChar_t *lastrow2
1062 ,AliHLTEtaRow etaclust,AliHLTTrackIndex *trackid
1066 // The method is a part of the fast hough transform.
1067 // It fills one row of the hough space.
1068 // It is called by FillCluster() method inside the
1069 // loop over alpha2 bins
1071 for(Int_t bin=binx1;bin<=binx2;bin++)
1073 if(ngaps2[bin] < MAX_N_GAPS) {
1074 if(i < lastrow2[bin] && i > currentrow2[bin]) {
1075 ngaps2[bin] += (i-currentrow2[bin]-1);
1079 if(i < lastrow2[bin] && i >= currentrow2[bin]) {
1080 for(UInt_t t=0;t<(MaxTrack-1); t++)
1082 Int_t label = etaclust.fMcLabels[t];
1083 if(label == 0) break;
1085 Int_t tempbin2 = (Int_t)(bin/2);
1086 for(c=0; c<(MaxTrack-1); c++)
1087 if(trackid[tempbin2].fLabel[c] == label || trackid[tempbin2].fNHits[c] == 0)
1089 trackid[tempbin2].fLabel[c] = label;
1090 if(trackid[tempbin2].fCurrentRow[c] != i) {
1091 trackid[tempbin2].fNHits[c]++;
1092 trackid[tempbin2].fCurrentRow[c] = i;
1102 inline void AliHLTTPCHoughTransformerRow::FillCluster(UChar_t i,Int_t etaindex,AliHLTEtaRow *etaclust,Int_t ilastpatch,Int_t firstbinx,Int_t lastbinx,Int_t nbinx,Int_t firstbiny)
1104 // The method is a part of the fast hough transform.
1105 // It fills a TPC cluster into the hough space.
1107 UChar_t *ngaps = fGapCount[etaindex];
1108 UChar_t *currentrow = fCurrentRowCount[etaindex];
1109 UChar_t *lastrow = fTrackLastRow;
1110 UChar_t *prevbin = fPrevBin[etaindex];
1111 UChar_t *nextbin = fNextBin[etaindex];
1112 UChar_t *nextrow = fNextRow[etaindex];
1114 AliHLTTrackIndex *trackid = fTrackID[etaindex];
1117 //Do the transformation:
1118 AliHLTPadHoughParams *startparams = &fStartPadParams[(Int_t)i][etaclust[etaindex].fStartPad];
1119 AliHLTPadHoughParams *endparams = &fEndPadParams[(Int_t)i][etaclust[etaindex].fEndPad];
1121 Float_t alpha1 = startparams->fAlpha;
1122 Float_t deltaalpha1 = startparams->fDeltaAlpha;
1123 Float_t alpha2 = endparams->fAlpha;
1124 Float_t deltaalpha2 = endparams->fDeltaAlpha;
1126 Int_t firstbin1 = startparams->fFirstBin;
1127 Int_t firstbin2 = endparams->fFirstBin;
1128 Int_t firstbin = firstbin1;
1129 if(firstbin>firstbin2) firstbin = firstbin2;
1131 Int_t lastbin1 = startparams->fLastBin;
1132 Int_t lastbin2 = endparams->fLastBin;
1133 Int_t lastbin = lastbin1;
1134 if(lastbin<lastbin2) lastbin = lastbin2;
1136 alpha1 += (firstbin-firstbiny)*deltaalpha1;
1137 alpha2 += (firstbin-firstbiny)*deltaalpha2;
1139 //Fill the histogram along the alpha2 range
1140 if(ilastpatch == -1) {
1141 for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
1143 Int_t binx1 = 1 + (Int_t)alpha1;
1144 if(binx1>lastbinx) continue;
1145 if(binx1<firstbinx) binx1 = firstbinx;
1146 Int_t binx2 = 1 + (Int_t)alpha2;
1147 if(binx2<firstbinx) continue;
1148 if(binx2>lastbinx) binx2 = lastbinx;
1151 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::TransformCircle()","")
1152 <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<etaclust[etaindex].fStartPad<<" "<<etaclust[etaindex].fEndPad<<ENDLOG;
1155 Int_t tempbin = b*nbinx;
1156 UChar_t *ngaps2 = ngaps + tempbin;
1157 UChar_t *currentrow2 = currentrow + tempbin;
1158 UChar_t *lastrow2 = lastrow + tempbin;
1160 Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2));
1161 AliHLTTrackIndex *trackid2 = trackid + tempbin2;
1163 FillClusterRow(i,binx1,binx2,ngaps2,currentrow2,lastrow2
1165 ,etaclust[etaindex],trackid2
1171 for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
1173 Int_t binx1 = 1 + (Int_t)alpha1;
1174 if(binx1>lastbinx) continue;
1175 if(binx1<firstbinx) binx1 = firstbinx;
1176 Int_t binx2 = 1 + (Int_t)alpha2;
1177 if(binx2<firstbinx) continue;
1178 if(binx2>lastbinx) binx2 = lastbinx;
1181 LOG(AliHLTTPCLog::kWarning,"AliHLTTPCHoughTransformerRow::TransformCircle()","")
1182 <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<etaclust[etaindex].fStartPad<<" "<<etaclust[etaindex].fEndPad<<ENDLOG;
1185 if(nextrow[b] > b) {
1186 Int_t deltab = (nextrow[b] - b - 1);
1188 alpha1 += deltaalpha1*deltab;
1189 alpha2 += deltaalpha2*deltab;
1192 Int_t tempbin = b*nbinx;
1193 binx1 = (UInt_t)nextbin[binx1 + tempbin];
1194 binx2 = (UInt_t)prevbin[binx2 + tempbin];
1195 if(binx2<binx1) continue;
1196 UChar_t *ngaps2 = ngaps + tempbin;
1197 UChar_t *currentrow2 = currentrow + tempbin;
1198 UChar_t *lastrow2 = lastrow + tempbin;
1200 Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2));
1201 AliHLTTrackIndex *trackid2 = trackid + tempbin2;
1203 FillClusterRow(i,binx1,binx2,ngaps2,currentrow2,lastrow2
1205 ,etaclust[etaindex],trackid2
1213 inline void AliHLTTPCHoughTransformerRow::FillClusterMCLabels(AliHLTTPCDigitData digpt,AliHLTEtaRow *etaclust)
1215 // The method is a part of the fast hough transform.
1216 // It fills the MC labels of a TPC cluster into a
1217 // special hough space array.
1218 for(Int_t t=0; t<3; t++)
1220 Int_t label = digpt.fTrackID[t];
1221 if(label < 0) break;
1223 for(c=0; c<(MaxTrack-1); c++)
1224 if(etaclust->fMcLabels[c] == label || etaclust->fMcLabels[c] == 0)
1227 etaclust->fMcLabels[c] = label;
1232 void AliHLTTPCHoughTransformerRow::SetTransformerArrays(AliHLTTPCHoughTransformerRow *tr)
1234 // In case of sequential filling of the hough space, the method is used to
1235 // transmit the pointers to the hough arrays from one transformer to the
1238 fGapCount = tr->fGapCount;
1239 fCurrentRowCount = tr->fCurrentRowCount;
1241 fTrackID = tr->fTrackID;
1243 fTrackNRows = tr->fTrackNRows;
1244 fTrackFirstRow = tr->fTrackFirstRow;
1245 fTrackLastRow = tr->fTrackLastRow;
1246 fInitialGapCount = tr->fInitialGapCount;
1248 fPrevBin = tr->fPrevBin;
1249 fNextBin = tr->fNextBin;
1250 fNextRow = tr->fNextRow;
1252 fStartPadParams = tr->fStartPadParams;
1253 fEndPadParams = tr->fEndPadParams;
1255 fLUTforwardZ = tr->fLUTforwardZ;
1256 fLUTbackwardZ = tr->fLUTbackwardZ;
1258 fParamSpace = tr->fParamSpace;