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