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