Bugfix in AliL3Hough::FindTrackCandidates; when track->SetEta, one has to
[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 eta_index,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   //Check if the track is leaving the sector at some point
97   Float_t maxrow=300;
98   Double_t angle=AliL3Transform::Pi()/18;
99   track->CalculateEdgePoint(angle);
100   if(!track->IsPoint())
101     {
102       track->CalculateEdgePoint(-1.*angle);
103       if(track->IsPoint())
104         maxrow = track->GetPointX();
105     }
106   else
107     maxrow = track->GetPointX();
108
109   for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
110     {
111       if(AliL3Transform::Row2X(padrow) > maxrow) break;//The track has left this slice
112       rows_crossed++;
113       Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
114       if(!track->GetCrossingPoint(padrow,xyz))  
115         {
116           continue;
117         }
118       
119       AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
120       AliL3Transform::Local2Raw(xyz,sector,row);
121       npixs=0;
122       
123       
124       //Get the timebins for this pad
125       AliL3DigitRowData *tempPt = fRowPointers[prow];
126       if(!tempPt) 
127         {
128           printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
129           continue;
130         }
131       
132       //Look at both sides of the pad:
133       for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
134         {
135           AliL3DigitData *digPt = tempPt->fDigitData;
136           for(UInt_t j=0; j<tempPt->fNDigit; j++)
137             {
138               UChar_t pad = digPt[j].fPad;
139               if(pad < p) continue;
140               if(pad > p) break;
141               UShort_t time = digPt[j].fTime;
142               Double_t eta = AliL3Transform::GetEta(padrow,pad,time);
143               Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);
144               if(pixel_index > eta_index) continue;
145               if(pixel_index != eta_index) break;
146               total_charge += digPt[j].fCharge;
147               if(remove)
148                 digPt[j].fCharge = 0; //Erase the track from image
149               npixs++;
150             }
151         }
152             
153       if(npixs > 1)//At least 2 digits on this padrow
154         {
155           nrow++;
156         }         
157     }
158   if(remove)
159     return kTRUE;
160   
161   if(nrow >= rows_crossed - fNumOfRowsToMiss)//this was a good track
162     {
163       if(fRemoveFoundTracks)
164         LookInsideRoad(track,eta_index,kTRUE);
165       return kTRUE;
166     }
167   else
168     return kFALSE;
169 }
170
171 void AliL3HoughEval::FindEta(AliL3TrackArray *tracks)
172 {
173   
174   Int_t sector,row;
175   Float_t xyz[3];
176   
177   Int_t ntracks = tracks->GetNTracks();
178   fEtaHistos = new AliL3Histogram1D*[ntracks];
179   
180   Char_t hname[100];
181   for(Int_t i=0; i<ntracks; i++)
182     {
183       sprintf(hname,"etahist_%d",i);
184       fEtaHistos[i] = new AliL3Histogram1D(hname,hname,100,0,1);
185     }
186   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
187   
188   for(Int_t ntr=0; ntr<ntracks; ntr++)
189     {
190       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(ntr);
191       if(!track) continue;
192       for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
193         {
194           Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
195           
196           if(!track->GetCrossingPoint(padrow,xyz))  
197             {
198               printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
199               continue;
200             }
201           
202           AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
203           AliL3Transform::Local2Raw(xyz,sector,row);
204           
205           //Get the timebins for this pad
206           AliL3DigitRowData *tempPt = fRowPointers[prow];
207           if(!tempPt) 
208             {
209               printf("AliL3HoughEval::LookInsideRoad : Zero data pointer\n");
210               continue;
211             }
212           
213           //Look at both sides of the pad:
214           for(Int_t p=(Int_t)rint(xyz[1])-fNumOfPadsToLook; p<=(Int_t)rint(xyz[1])+fNumOfPadsToLook; p++)
215             {
216               AliL3DigitData *digPt = tempPt->fDigitData;
217               for(UInt_t j=0; j<tempPt->fNDigit; j++)
218                 {
219                   UChar_t pad = digPt[j].fPad;
220                   
221                   if(pad < p) continue;
222                   if(pad > p) break;
223                   UShort_t time = digPt[j].fTime;
224                   Double_t eta = AliL3Transform::GetEta(padrow,pad,time);
225                   Int_t pixel_index = (Int_t)(eta/etaslice);
226                   if(pixel_index > track->GetEtaIndex()+1) continue;
227                   if(pixel_index < track->GetEtaIndex()-1) break;
228                   fEtaHistos[ntr]->Fill(eta,digPt[j].fCharge);
229                 }
230             }
231         }
232     }
233   
234   for(Int_t i=0; i<ntracks; i++)
235     {
236       AliL3Histogram1D *hist = fEtaHistos[i];
237       Int_t max_bin = hist->GetMaximumBin();
238       Double_t max_value = hist->GetBinContent(max_bin);
239       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
240       if(!track) continue;
241       if(hist->GetBinContent(max_bin-1)<max_value && hist->GetBinContent(max_bin+1)<max_value)
242         {
243           track->SetWeight((Int_t)max_value,kTRUE); 
244           track->SetEta(hist->GetBinCenter(max_bin));
245           track->SetNHits(track->GetWeight());
246         }
247       else
248         {
249           track->SetWeight(0);
250           tracks->Remove(i); //remove this track, because it was not a peak
251         }    
252     }
253   tracks->Compress();
254   
255   //for(Int_t i=0; i<ntracks; i++)
256   //delete fEtaHistos[i];
257   //delete []¬†fEtaHistos;
258 }
259
260 void AliL3HoughEval::DisplayEtaSlice(Int_t eta_index,AliL3Histogram *hist)
261 {
262   //Display the current raw data inside the (slice,patch)
263
264   if(!hist)
265     {
266       printf("AliL3HoughEval::DisplayEtaSlice : No input histogram!\n");
267       return;
268     }
269   
270   for(Int_t padrow = AliL3Transform::GetFirstRow(fPatch); padrow <= AliL3Transform::GetLastRow(fPatch); padrow++)
271     {
272       Int_t prow = padrow - AliL3Transform::GetFirstRow(fPatch);
273                   
274       AliL3DigitRowData *tempPt = fRowPointers[prow];
275       if(!tempPt) 
276         {
277           printf("AliL3HoughEval::DisplayEtaSlice : Zero data pointer\n");
278           continue;
279         }
280       
281       AliL3DigitData *digPt = tempPt->fDigitData;
282       if((Int_t)tempPt->fRow != padrow)
283         {
284           printf("\nAliL3HoughEval::DisplayEtaSlice : Mismatching padrows!!!\n");
285           return;
286         }
287       for(UInt_t j=0; j<tempPt->fNDigit; j++)
288         {
289           UChar_t pad = digPt[j].fPad;
290           UChar_t charge = digPt[j].fCharge;
291           UShort_t time = digPt[j].fTime;
292           if((Int_t)charge < fHoughTransformer->GetLowerThreshold() || (Int_t)charge > fHoughTransformer->GetUpperThreshold()) continue;
293           Float_t xyz[3];
294           Int_t sector,row;
295           AliL3Transform::Slice2Sector(fSlice,padrow,sector,row);
296           AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
297           Double_t eta = AliL3Transform::GetEta(xyz);
298           Int_t pixel_index = fHoughTransformer->GetEtaIndex(eta);//(Int_t)(eta/etaslice);
299           if(pixel_index != eta_index) continue;
300           hist->Fill(xyz[0],xyz[1],charge);
301         }
302     }
303   
304 }
305
306 #ifdef use_root
307 void AliL3HoughEval::CompareMC(AliL3TrackArray *tracks,Char_t *trackfile,Int_t threshold)
308 {
309   /*  
310   struct GoodTrack goodtracks[15000];
311   Int_t nt=0;
312   ifstream in(trackfile);
313   if(in)
314     {
315       printf("Reading good tracks from file %s\n",trackfile);
316       while (in>>goodtracks[nt].label>>goodtracks[nt].code>>
317              goodtracks[nt].px>>goodtracks[nt].py>>goodtracks[nt].pz>>
318              goodtracks[nt].pt>>goodtracks[nt].eta>>goodtracks[nt].nhits) 
319         {
320           nt++;
321           if (nt==15000) 
322             {
323               cerr<<"Too many good tracks"<<endl;
324               break;
325             }
326         }
327       if (!in.eof())
328         {
329           LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input file")
330             <<"Error in file reading"<<ENDLOG;
331           return;
332         }
333     }
334   else
335     {
336       LOG(AliL3Log::kError,"AliL3HoughEval::CompareMC","Input")
337         <<"No input trackfile "<<trackfile<<ENDLOG;
338     }
339   
340   Int_t *particles = new Int_t[fNEtaSegments];
341   Int_t *ftracks = new Int_t[fNEtaSegments];
342   for(Int_t i=0; i<fNEtaSegments; i++)
343     {
344       particles[i]=0;
345       ftracks[i]=0;
346     }
347   
348   TH1F *ptgood = new TH1F("ptgood","ptgood",5,0,2);
349   TH1F *ptfound = new TH1F("ptfound","ptgood",5,0,2);
350   TH1F *pteff = new TH1F("pteff","pteff",5,0,2);
351   TH1F *etafound = new TH1F("etafound","etafound",5,0,1);
352   TH1F *etagood = new TH1F("etagood","etagood",5,0,1);
353   TH1F *etaeff = new TH1F("etaeff","etaeff",5,0,1);
354   
355   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
356   for(Int_t i=0; i<tracks->GetNTracks(); i++)
357     {
358       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
359       if(!tr) continue;
360       if(tr->GetWeight()<threshold) continue;
361       Int_t trackindex = tr->GetEtaIndex();
362       if(trackindex <0 || trackindex >= fNEtaSegments) continue;
363       ftracks[trackindex]++;
364       ptfound->Fill(tr->GetPt());
365       etafound->Fill(tr->GetEta());
366     }
367   for(Int_t i=0; i<nt; i++)
368     {
369       if(goodtracks[i].nhits < 174) continue;
370       if(goodtracks[i].pt < 0.2) continue;
371       Int_t particleindex = (Int_t)(goodtracks[i].eta/etaslice);
372       if(particleindex < 0 || particleindex >= fNEtaSegments) continue;
373       particles[particleindex]++;
374       ptgood->Fill(goodtracks[i].pt);
375       etagood->Fill(goodtracks[i].eta);
376     }
377   
378   Double_t found=0;
379   Double_t good =0;
380   for(Int_t i=0; i<fNEtaSegments; i++)
381     {
382       //printf("Slice %d : Found tracks %d, good tracks %d\n",i,ftracks[i],particles[i]);
383       found += ftracks[i];
384       good += particles[i];
385     }
386   printf("And the total efficiency was: %f\n",found/good);
387
388   ptgood->Sumw2(); ptfound->Sumw2();
389   etagood->Sumw2(); etafound->Sumw2();
390   pteff->Divide(ptfound,ptgood,1,1,"b");
391   etaeff->Divide(etafound,etagood,1,1,"b");
392   TFile *file = TFile::Open("eff.root","RECREATE");
393   ptgood->Write();
394   ptfound->Write();
395   pteff->Write();
396   etafound->Write();
397   etagood->Write();
398   etaeff->Write();
399   file->Close();
400   
401   delete [] particles;
402   delete [] ftracks;
403   */  
404 }
405
406 #endif