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