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