d7d168ff2fd8756321217e299648f627f4a69987
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughEval.cxx
1 //$Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ASV 
5
6 #include <math.h>
7 #include <iostream.h>
8 #ifdef use_root
9 #include <TH1.h>
10 #include <TFile.h>
11 #endif
12 #include "AliL3MemHandler.h"
13 #include "AliL3TrackArray.h"
14 #include "AliL3Logging.h"
15 #include "AliL3HoughEval.h"
16 #include "AliL3HoughBaseTransformer.h"
17 #include "AliL3DigitData.h"
18 #include "AliL3HoughTrack.h"
19 #include "AliL3Transform.h"
20 #include "AliL3Histogram.h"
21 #include "AliL3Histogram1D.h"
22
23 //_____________________________________________________________
24 // AliL3HoughEval
25 //
26 // Evaluation class for tracklets produced by the Hough transform.
27
28 ClassImp(AliL3HoughEval)
29
30 AliL3HoughEval::AliL3HoughEval()
31 {
32   
33   fRemoveFoundTracks = kFALSE;
34   fNumOfPadsToLook = 1;
35   fNumOfRowsToMiss = 1;
36   fEtaHistos=0;
37   fRowPointers = 0;
38 }
39
40
41 AliL3HoughEval::~AliL3HoughEval()
42 {
43   fHoughTransformer = 0;
44   if(fRowPointers)
45     {
46       for(Int_t i=0; i<fNrows; i++)
47         fRowPointers[i] = 0;
48       delete [] fRowPointers;
49     }
50 }
51
52 void AliL3HoughEval::InitTransformer(AliL3HoughBaseTransformer *transformer)
53 {
54   fHoughTransformer = transformer;
55   fSlice = fHoughTransformer->GetSlice();
56   fPatch = fHoughTransformer->GetPatch();
57   fNrows = AliL3Transform::GetLastRow(fPatch) - AliL3Transform::GetFirstRow(fPatch) + 1;
58   fNEtaSegments = fHoughTransformer->GetNEtaSegments();
59   fEtaMin = fHoughTransformer->GetEtaMin();
60   fEtaMax = fHoughTransformer->GetEtaMax();
61   GenerateLUT();
62 }
63
64 void AliL3HoughEval::GenerateLUT()
65 {
66   //Generate a Look-up table, to limit the access to raw data
67   
68   if(!fRowPointers)
69     fRowPointers = new AliL3DigitRowData*[fNrows];
70
71   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer->GetDataPointer();
72   if(!tempPt)
73     printf("\nAliL3HoughEval::GenerateLUT : Zero data pointer\n");
74   
75   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
76     {
77       Int_t prow = i - AliL3Transform::GetFirstRow(fPatch);
78       fRowPointers[prow] = tempPt;
79       AliL3MemHandler::UpdateRowPointer(tempPt);
80     }
81   
82 }
83
84 Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t &nrows_crossed,Bool_t remove)
85 {
86   //Look at rawdata along the road specified by the track candidates.
87   //If track is good, return true, if not return false.
88   
89   Int_t sector,row;
90   
91   Int_t nrow=0,npixs=0;//,rows_crossed=0;
92   Float_t xyz[3];
93   
94   Int_t total_charge=0;//total charge along the road
95     
96   for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
97     {
98       Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
99       if(!track->GetCrossingPoint(padrow,xyz))  
100         {
101           continue;
102         }
103       
104       AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
105       AliL3Transform::Local2Raw(xyz,sector,row);
106
107       npixs=0;
108       
109       //Get the timebins for this pad
110       AliL3DigitRowData *tempPt = fRowPointers[prow];
111       if(!tempPt) 
112         {
113           printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
114           continue;
115         }
116       
117       //Look at both sides of the pad:
118       for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
119         {
120           AliL3DigitData *digPt = tempPt->fDigitData;
121           for(UInt_t j=0; j<tempPt->fNDigit; j++)
122             {
123               UChar_t pad = digPt[j].fPad;
124               Int_t charge = digPt[j].fCharge;
125               if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
126               if(pad < p) continue;
127               if(pad > p) break;
128               UShort_t time = digPt[j].fTime;
129               Double_t eta = AliL3Transform::GetEta(padrow,pad,time);
130               Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);
131               if(pixel_index != track->GetEtaIndex()) continue;
132               total_charge += digPt[j].fCharge;
133               if(remove)
134                 digPt[j].fCharge = 0; //Erase the track from image
135               npixs++;
136             }
137         }
138             
139       if(npixs > 1)//At least 2 digits on this padrow
140         {
141           nrow++;
142         }         
143     }
144   if(remove)
145     return kTRUE;
146   
147   nrows_crossed += nrow; //Update the number of rows crossed.
148   
149   if(nrow >= AliL3Transform::GetNRows(fPatch) - fNumOfRowsToMiss)//this was a good track
150     {
151       if(fRemoveFoundTracks)
152         {
153           Int_t dummy=0;
154           LookInsideRoad(track,dummy,kTRUE);
155         }
156       return kTRUE;
157     }
158   else
159     return kFALSE;
160 }
161
162 void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
163 {
164   
165   Int_t sector,row;
166   Float_t xyz[3];
167   
168   Int_t ntracks = tracks->GetNTracks();
169   fEtaHistos = new AliL3Histogram1D*[ntracks];
170   
171   Char_t hname[100];
172   for(Int_t i=0; i<ntracks; i++)
173     {
174       sprintf(hname,"etahist_%d",i);
175       fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
176     }
177   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
178   
179   for(Int_t ntr=0; ntr<ntracks; ntr++)
180     {
181       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
182       if(!track) continue;
183       for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
184         {
185           Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
186           
187           if(!track->GetCrossingPoint(padrow,xyz))  
188             {
189               printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
190               continue;
191             }
192           
193           AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
194           AliL3Transform::Local2Raw(xyz,sector,row);
195           
196           //Get the timebins for this pad
197           AliL3DigitRowData *tempPt = fRowPointers[prow];
198           if(!tempPt) 
199             {
200               printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
201               continue;
202             }
203           
204           //Look at both sides of the pad:
205           for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
206             {
207               AliL3DigitData *digPt = tempPt->fDigitData;
208               for(UInt_t j=0; j<tempPt->fNDigit; j++)
209                 {
210                   UChar_t pad = digPt[j].fPad;
211                   Int_t charge = digPt[j].fCharge;
212                   if(charge <= fHoughTransformer->GetLowerThreshold()) continue;
213                   if(pad < p) continue;
214                   if(pad > p) break;
215                   UShort_t time = digPt[j].fTime;
216                   Double_t eta = AliL3Transform::GetEta(padrow,pad,time);
217                   Int_t pixel_index = (Int_t)(eta/etaslice);
218                   if(pixel_index > track->GetEtaIndex()+1) continue;
219                   if(pixel_index < track->GetEtaIndex()-1) break;
220                   fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
221                 }
222             }
223         }
224     }
225   
226   for(Int_t i=0; i<ntracks; i++)
227     {
228       AliL3Histogram1D *hist = fEtaHistos[i];
229       Int_t max_bin = hist->GetMaximumBin();
230       Double_t max_value = hist->GetBinContent(max_bin);
231       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
232       if(!track) continue;
233       if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
234         {
235           track->SetWeight((Int_t)max_value,kTRUE); 
236           track->SetEta(hist->GetBinCenter(max_bin));
237           track->SetNHits(track->GetWeight());
238         }
239       else
240         {
241           track->SetWeight(0);
242           tracks->Remove(i); //remove this track, because it was not a peak
243         }    
244     }
245   tracks->Compress();
246   
247   //for(Int_t i=0; i<ntracks; i++)
248   //delete fEtaHistos[i];
249   //delete []¬†fEtaHistos;
250 }
251
252 void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
253 {
254   //Display the current raw data inside the (slice,patch)
255
256   if(!hist)
257     {
258       printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
259       return;
260     }
261   
262   for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
263     {
264       Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
265                   
266       AliL3DigitRowData *tempPt = fRowPointers[prow];
267       if(!tempPt) 
268         {
269           printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\n");
270           continue;
271         }
272       
273       AliL3DigitData *digPt = tempPt->fDigitData;
274       if((Int_t)tempPt->fRow != padrow)
275         {
276           printf("\nAliL3HoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
277           return;
278         }
279       for(UInt_t j=0; j<tempPt->fNDigit; j++)
280         {
281           UChar_t pad = digPt[j].fPad;
282           UChar_t charge = digPt[j].fCharge;
283           UShort_t time = digPt[j].fTime;
284           if((Int_t)charge <= fHoughTransformer->GetLowerThreshold() || (Int_t)charge >= fHoughTransformer->GetUpperThreshold()) continue;
285           Float_t xyz[3];
286           Int_t sector,row;
287           AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
288           AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
289           Double_t eta = AliL3Transform::GetEta(xyz);
290           Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
291           if(pixel_index != eta_index) continue;
292           hist->Fill(xyz[0],xyz[1],charge);
293         }
294     }
295   
296 }
297
298 #ifdef use_root
299 void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile,Int_t threshold)
300 {
301   /*  
302   struct GoodTrack goodtracks[15000];
303   Int_t nt=0;
304   ifstream in(trackfile);
305   if(in)
306     {
307       printf("Reading good tracks from file %s\n",trackfile);
308       while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
309              goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
310              goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits) 
311         {
312           nt++;
313           if (nt==15000) 
314             {
315               cerr<<"Too many good tracks"<<endl;
316               break;
317             }
318         }
319       if (!in.eof())
320         {
321           LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input file")
322             <<"Error in file reading"<<ENDLOG;
323           return;
324         }
325     }
326   else
327     {
328       LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input")
329         <<"No input trackfile "<<trackfile<<ENDLOG;
330     }
331   
332   Int_t *particles = new Int_t[fNEtaSegments];
333   Int_t *ftracks = new Int_t[fNEtaSegments];
334   for(Int_t i=0; i<fNEtaSegments; i++)
335     {
336       particles[i]=0;
337       ftracks[i]=0;
338     }
339   
340   TH1F *ptgood = new TH1F("ptgood","ptgood",5,0,2);
341   TH1F *ptfound = new TH1F("ptfound","ptgood",5,0,2);
342   TH1F *pteff = new TH1F("pteff","pteff",5,0,2);
343   TH1F *etafound = new TH1F("etafound","etafound",5,0,1);
344   TH1F *etagood = new TH1F("etagood","etagood",5,0,1);
345   TH1F *etaeff = new TH1F("etaeff","etaeff",5,0,1);
346   
347   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
348   for(Int_t i=0; i<tracks->GetNTracks(); i++)
349     {
350       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
351       if(!tr) continue;
352       if(tr->GetWeight()<threshold) continue;
353       Int_t trackindex = tr->GetEtaIndex();
354       if(trackindex <0 || trackindex >= fNEtaSegments) continue;
355       ftracks[trackindex]++;
356       ptfound->Fill(tr->GetPt());
357       etafound->Fill(tr->GetEta());
358     }
359   for(Int_t i=0; i<nt; i++)
360     {
361       if(goodtracks[i].nhits < 174) continue;
362       if(goodtracks[i].pt < 0.2) continue;
363       Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
364       if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
365       particles[particleindex]++;
366       ptgood->Fill(goodtracks[i].pt);
367       etagood->Fill(goodtracks[i].eta);
368     }
369   
370   Double_t found=0;
371   Double_t good =0;
372   for(Int_t i=0; i<fNEtaSegments; i++)
373     {
374       //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
375       found += ftracks[i];
376       good += particles[i];
377     }
378   printf("And the total efficiency was: %f\n",found/good);
379
380   ptgood->Sumw2(); ptfound->Sumw2();
381   etagood->Sumw2(); etafound->Sumw2();
382   pteff->Divide(ptfound,ptgood,1,1,"b");
383   etaeff->Divide(etafound,etagood,1,1,"b");
384   TFile *file = TFile::Open("eff.root","RECREATE");
385   ptgood->Write();
386   ptfound->Write();
387   pteff->Write();
388   etafound->Write();
389   etagood->Write();
390   etaeff->Write();
391   file->Close();
392   
393   delete [] particles;
394   delete [] ftracks;
395   */  
396 }
397
398 #endif