3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ASV
6 #include "AliL3StandardIncludes.h"
9 #include "AliL3Logging.h"
10 #include "AliL3HoughMerger.h"
11 #include "AliL3HoughIntMerger.h"
12 #include "AliL3HoughGlobalMerger.h"
13 #include "AliL3Histogram.h"
14 #include "AliL3Hough.h"
15 #include "AliL3HoughTransformer.h"
16 #include "AliL3HoughClusterTransformer.h"
17 #include "AliL3HoughTransformerLUT.h"
18 #include "AliL3HoughTransformerVhdl.h"
19 #include "AliL3HoughMaxFinder.h"
21 #include "AliL3FileHandler.h"
23 #include "AliL3MemHandler.h"
25 #include "AliL3DataHandler.h"
26 #include "AliL3DigitData.h"
27 #include "AliL3HoughEval.h"
28 #include "AliL3Transform.h"
29 #include "AliL3TrackArray.h"
30 #include "AliL3HoughTrack.h"
38 //_____________________________________________________________
41 // Interface class for the Hough transform
43 // Example how to use:
45 // AliL3Hough *hough = new AliL3Hough(path,kTRUE,NumberOfEtaSegments);
46 // hough->ReadData(slice);
47 // hough->Transform();
48 // hough->FindTrackCandidates();
50 // AliL3TrackArray *tracks = hough->GetTracks(patch);
56 AliL3Hough::AliL3Hough()
61 fAddHistograms = kFALSE;
62 fDoIterative = kFALSE;
63 fWriteDigits = kFALSE;
67 fHoughTransformer = 0;
80 SetTransformerParams();
85 AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv)
91 fNEtaSegments = n_eta_segments;
92 fAddHistograms = kFALSE;
93 fDoIterative = kFALSE;
94 fWriteDigits = kFALSE;
99 AliL3Hough::~AliL3Hough()
111 delete fGlobalMerger;
114 void AliL3Hough::CleanUp()
118 for(Int_t i=0; i<fNPatches; i++)
120 if(fTracks[i]) delete fTracks[i];
121 if(fEval[i]) delete fEval[i];
122 if(fHoughTransformer[i]) delete fHoughTransformer[i];
123 if(fMemHandler[i]) delete fMemHandler[i];
127 if(fTracks) delete [] fTracks;
128 if(fEval) delete [] fEval;
129 if(fHoughTransformer) delete [] fHoughTransformer;
130 if(fMemHandler) delete [] fMemHandler;
134 void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv)
138 fNEtaSegments = n_eta_segments;
139 fWriteDigits = kFALSE;
143 Init(); //do the rest
146 void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
149 fAddHistograms = addhists;
151 AliL3Transform::Init(fPath,!fBinary);
152 fNPatches = AliL3Transform::GetNPatches();
154 fHoughTransformer = new AliL3HoughBaseTransformer*[fNPatches];
155 fMemHandler = new AliL3MemHandler*[fNPatches];
156 fTracks = new AliL3TrackArray*[fNPatches];
157 fEval = new AliL3HoughEval*[fNPatches];
159 for(Int_t i=0; i<fNPatches; i++)
161 switch (fVersion){ //choose Transformer
163 fHoughTransformer[i] = new AliL3HoughTransformerLUT(0,i,fNEtaSegments);
166 fHoughTransformer[i] = new AliL3HoughClusterTransformer(0,i,fNEtaSegments);
169 fHoughTransformer[i] = new AliL3HoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
172 fHoughTransformer[i] = new AliL3HoughTransformer(0,i,fNEtaSegments);
175 fHoughTransformer[i]->CreateHistograms(fNBinX,fLowPt,fNBinY,-fPhi,fPhi);
176 fHoughTransformer[i]->SetLowerThreshold(fThreshold);
178 LOG(AliL3Log::kInformational,"AliL3Hough::Init","Version")
179 <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
181 fEval[i] = new AliL3HoughEval();
182 fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
184 fMemHandler[i] = new AliL3DataHandler();
188 fMemHandler[i] = new AliL3FileHandler();
191 Char_t filename[1024];
192 sprintf(filename,"%s/digitfile.root",fPath);
193 fMemHandler[i]->SetAliInput(filename);
197 fMemHandler[i] = new AliL3MemHandler();
201 fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",1000);
202 fMerger = new AliL3HoughMerger(fNPatches);
203 fInterMerger = new AliL3HoughIntMerger();
207 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
209 //Process all slices [minslice,maxslice].
210 fGlobalMerger = new AliL3HoughGlobalMerger(minslice,maxslice);
212 for(Int_t i=minslice; i<=maxslice; i++)
218 FindTrackCandidates();
220 fGlobalMerger->FillTracks(fTracks[0],i);
224 void AliL3Hough::ReadData(Int_t slice,Int_t eventnr)
226 //Read data from files, binary or root.
228 fCurrentSlice = slice;
229 for(Int_t i=0; i<fNPatches; i++)
231 fMemHandler[i]->Free();
233 AliL3DigitRowData *digits =0;
235 fMemHandler[i]->Init(slice,i);
236 if(fBinary)//take input data from binary files
239 sprintf(name,"%sdigits_c8_%d_%d.raw",fPath,slice,i);
241 sprintf(name,"%sdigits_%d_%d.raw",fPath,slice,i);
243 fMemHandler[i]->SetBinaryInput(name);
244 digits = (AliL3DigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
245 fMemHandler[i]->CloseBinaryInput();
247 else //read data from root file
250 digits=(AliL3DigitRowData *)fMemHandler[i]->AliDigits2Memory(ndigits,eventnr);
251 fMemHandler[i]->FreeDigitsTree();
253 cerr<<"You cannot read from rootfile now"<<endl;
257 //set input data and init transformer
258 fHoughTransformer[i]->SetInputData(ndigits,digits);
259 fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
263 void AliL3Hough::Transform(Int_t row_range)
265 //Transform all data given to the transformer within the given slice
266 //(after ReadData(slice))
268 Double_t initTime,cpuTime;
269 initTime = GetCpuTime();
270 for(Int_t i=0; i<fNPatches; i++)
272 fHoughTransformer[i]->Reset();//Reset the histograms
274 fHoughTransformer[i]->TransformCircle();
276 fHoughTransformer[i]->TransformCircleC(row_range);
278 cpuTime = GetCpuTime() - initTime;
279 LOG(AliL3Log::kInformational,"AliL3Hough::Transform()","Timing")
280 <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
283 void AliL3Hough::MergePatches()
285 if(fAddHistograms) //Nothing to merge here
287 fMerger->MergePatches(kTRUE);
290 void AliL3Hough::MergeInternally()
293 fInterMerger->FillTracks(fTracks[0]);
295 fInterMerger->FillTracks(fMerger->GetOutTracks());
297 fInterMerger->MMerge();
300 void AliL3Hough::ProcessSliceIter()
302 //Process current slice (after ReadData(slice)) iteratively.
306 for(Int_t i=0; i<fNPatches; i++)
309 fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
314 for(Int_t i=0; i<10; i++)
319 AliL3HoughBaseTransformer *tr = fHoughTransformer[0];
320 for(Int_t j=0; j<fNEtaSegments; j++)
322 AliL3Histogram *hist = tr->GetHistogram(j);
323 if(hist->GetNEntries()==0) continue;
324 fPeakFinder->Reset();
325 fPeakFinder->SetHistogram(hist);
326 fPeakFinder->FindAbsMaxima();
327 AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[0]->NextTrack();
328 track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
329 track->SetEtaIndex(j);
330 track->SetEta(tr->GetEta(j,fCurrentSlice));
331 for(Int_t k=0; k<fNPatches; k++)
333 fEval[i]->SetNumOfPadsToLook(2);
334 fEval[i]->SetNumOfRowsToMiss(2);
335 fEval[i]->RemoveFoundTracks();
337 if(!fEval[i]->LookInsideRoad(track,nrows))
339 fTracks[0]->Remove(fTracks[0]->GetNTracks()-1);
340 fTracks[0]->Compress();
350 void AliL3Hough::ProcessPatchIter(Int_t patch)
352 //Process patch in a iterative way.
353 //transform + peakfinding + evaluation + transform +...
355 Int_t num_of_tries = 5;
356 AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
357 AliL3TrackArray *tracks = fTracks[patch];
359 AliL3HoughEval *ev = fEval[patch];
360 ev->InitTransformer(tr);
361 //ev->RemoveFoundTracks();
362 ev->SetNumOfRowsToMiss(3);
363 ev->SetNumOfPadsToLook(2);
364 AliL3Histogram *hist;
365 for(Int_t t=0; t<num_of_tries; t++)
368 tr->TransformCircle();
369 for(Int_t i=0; i<fNEtaSegments; i++)
371 hist = tr->GetHistogram(i);
372 if(hist->GetNEntries()==0) continue;
373 fPeakFinder->Reset();
374 fPeakFinder->SetHistogram(hist);
375 fPeakFinder->FindAbsMaxima();
376 //fPeakFinder->FindPeak1();
377 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
378 track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
379 track->SetEtaIndex(i);
380 track->SetEta(tr->GetEta(i,fCurrentSlice));
382 if(!ev->LookInsideRoad(track,nrows))
384 tracks->Remove(tracks->GetNTracks()-1);
390 LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
391 <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
394 void AliL3Hough::AddAllHistograms()
396 //Add the histograms within one etaslice.
397 //Resulting histogram are in patch=0.
399 Double_t initTime,cpuTime;
400 initTime = GetCpuTime();
401 for(Int_t i=0; i<fNEtaSegments; i++)
403 AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
404 for(Int_t j=1; j<fNPatches; j++)
406 AliL3Histogram *hist = fHoughTransformer[j]->GetHistogram(i);
410 fAddHistograms = kTRUE;
411 cpuTime = GetCpuTime() - initTime;
412 LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistograms()","Timing")
413 <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
416 void AliL3Hough::FindTrackCandidates()
418 //Look for peaks in histograms, and find the track candidates
422 n_patches = 1; //Histograms have been added.
424 n_patches = fNPatches;
426 Double_t initTime,cpuTime;
427 initTime = GetCpuTime();
429 for(Int_t i=0; i<n_patches; i++)
431 AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
434 for(Int_t j=0; j<fNEtaSegments; j++)
436 AliL3Histogram *hist = tr->GetHistogram(j);
437 if(hist->GetNEntries()==0) continue;
438 fPeakFinder->Reset();
439 fPeakFinder->SetHistogram(hist);
440 //fPeakFinder->FindPeak1(3,1);
441 fPeakFinder->FindMaxima(0,0); //Simple maxima finder
442 //fPeakFinder->FindAbsMaxima();
443 for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
445 if(fPeakFinder->GetWeight(k) == 0) continue;
446 AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
447 track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
448 track->SetEtaIndex(j);
449 track->SetEta(tr->GetEta(j,fCurrentSlice));
450 track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
455 cpuTime = GetCpuTime() - initTime;
456 LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
457 <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
460 void AliL3Hough::InitEvaluate()
462 //Pass the transformer objects to the AliL3HoughEval objects:
463 //This will provide the evaluation objects with all the necessary
464 //data and parameters it needs.
466 for(Int_t i=0; i<fNPatches; i++)
467 fEval[i]->InitTransformer(fHoughTransformer[i]);
470 Int_t AliL3Hough::Evaluate(Int_t road_width,Int_t nrowstomiss)
472 //Evaluate the tracks, by looking along the road in the raw data.
473 //If track does not cross all padrows - rows2miss, it is removed from the arrray.
474 //If histograms were not added, the check is done locally in patch,
475 //meaning that nrowstomiss is the number of padrows the road can miss with respect
476 //to the number of rows in the patch.
477 //If the histograms were added, the comparison is done globally in the _slice_,
478 //meaing that nrowstomiss is the number of padrows the road can miss with
479 //respect to the total number of padrows in the slice.
481 //Return value = number of tracks which were removed (only in case of fAddHistograms)
485 LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
486 <<"No tracks to work with..."<<ENDLOG;
490 Int_t removed_tracks=0;
491 AliL3TrackArray *tracks=0;
496 for(Int_t i=0; i<tracks->GetNTracks(); i++)
498 AliL3Track *track = tracks->GetCheckedTrack(i);
504 for(Int_t i=0; i<fNPatches; i++)
505 EvaluatePatch(i,road_width,nrowstomiss);
507 //Here we check the tracks globally;
508 //how many good rows (padrows with signal)
509 //did it cross in the slice
513 for(Int_t j=0; j<tracks->GetNTracks(); j++)
515 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
517 if(track->GetNHits() < AliL3Transform::GetNRows() - nrowstomiss)
527 return removed_tracks;
530 void AliL3Hough::EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss)
534 fEval[i]->InitTransformer(fHoughTransformer[i]);
535 fEval[i]->SetNumOfPadsToLook(road_width);
536 fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
537 //fEval[i]->RemoveFoundTracks();
539 AliL3TrackArray *tracks=0;
547 for(Int_t j=0; j<tracks->GetNTracks(); j++)
549 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
552 LOG(AliL3Log::kWarning,"AliL3Hough::EvaluatePatch","Track array")
553 <<"Track object missing!"<<ENDLOG;
557 Bool_t result = fEval[i]->LookInsideRoad(track,nrows);
560 Int_t pre=track->GetNHits();
561 track->SetNHits(pre+nrows);
563 //else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
576 fMerger->FillTracks(tracks,i); //Copy tracks to the track merger
582 void AliL3Hough::MergeEtaSlices()
584 //Merge tracks found in neighbouring eta slices.
585 //Removes the track with the lower weight.
587 AliL3TrackArray *tracks = fTracks[0];
590 cerr<<"AliL3Hough::MergeEtaSlices : No tracks "<<endl;
593 for(Int_t j=0; j<tracks->GetNTracks(); j++)
595 AliL3HoughTrack *track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
596 if(!track1) continue;
597 for(Int_t k=j+1; k<tracks->GetNTracks(); k++)
599 AliL3HoughTrack *track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(k);
600 if(!track2) continue;
601 if(abs(track1->GetEtaIndex() - track2->GetEtaIndex()) != 1) continue;
602 if(track1->GetKappa() == track2->GetKappa() && track1->GetPsi() == track2->GetPsi())
604 cout<<"Merging track in slices "<<track1->GetEtaIndex()<<" "<<track2->GetEtaIndex()<<endl;
605 if(track1->GetWeight() > track2->GetWeight())
615 void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
617 //Write the tracks in slice
619 AliL3MemHandler *mem = new AliL3MemHandler();
623 sprintf(fname,"%s/tracks_ho_%d.raw",path,slice);
624 mem->SetBinaryOutput(fname);
625 mem->TrackArray2Binary(fTracks[0]);
626 mem->CloseBinaryOutput();
630 for(Int_t i=0; i<fNPatches; i++)
632 sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,slice,i);
633 mem->SetBinaryOutput(fname);
634 mem->TrackArray2Binary(fTracks[i]);
635 mem->CloseBinaryOutput();
642 void AliL3Hough::WriteDigits(Char_t *outfile)
645 //Write the current data to a new rootfile.
647 for(Int_t i=0; i<fNPatches; i++)
649 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer[i]->GetDataPointer();
650 fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
653 cerr<<"AliL3Hough::WriteDigits : You need to compile with AliROOT!"<<endl;
658 Double_t AliL3Hough::GetCpuTime()
660 //Return the Cputime in seconds.
662 gettimeofday( &tv, NULL );
663 return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
664 //return (Double_t)(clock()) / CLOCKS_PER_SEC;