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