18fa208aed6bff4a54201b4019dd0a5163087272
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughEval.cxx
1 //Author:        Anders Strand Vestbo
2 //Last Modified: 28.6.01
3
4 #include <TClonesArray.h>
5 #include <TH2.h>
6 #include <TH1.h>
7 #include <TFile.h>
8 #include <AliRun.h>
9 #include <TParticle.h>
10 #include <TTree.h>
11
12 #include "AliTPCParam.h"
13 #include "AliSimDigits.h"
14 #include "AliL3TrackArray.h"
15 #include "AliL3Transform.h"
16 #include "AliL3HoughTransformer.h"
17 #include "AliL3Defs.h"
18 #include "AliL3HoughTrack.h"
19 #include "AliL3HoughEval.h"
20
21 ClassImp(AliL3HoughEval)
22
23
24 AliL3HoughEval::AliL3HoughEval()
25 {
26   //Default constructor
27   fTransform = new AliL3Transform();
28   fHoughTransformer = 0;
29   fNumOfRowsToMiss = 1;
30   fNumOfPadsToLook = 1;
31 }
32
33
34 AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
35 {
36   //Constructor
37   fHoughTransformer = transformer;
38   fTransform = new AliL3Transform();
39   fNumOfRowsToMiss = 1;
40   fNumOfPadsToLook = 1;
41 }
42
43
44 AliL3HoughEval::~AliL3HoughEval()
45 {
46   //Destructor
47   if(fTransform)
48     delete fTransform;
49   if(fMcTrackTable)
50     delete [] fMcTrackTable;
51 }
52
53
54 Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
55 {
56   //Look at rawdata along the road specified by the track candidates.
57   
58   if(!fHoughTransformer)
59     {
60       printf("\nAliL3HoughEval: No transformer object\n");
61       return kFALSE;
62     }
63   
64   Int_t patch = fHoughTransformer->fPatch;
65   Int_t slice = fHoughTransformer->fSlice;
66
67   
68   Int_t sector,row;
69   Int_t lut_index;
70   Char_t **track_lut = fHoughTransformer->fTrackTable;
71   Int_t nrow=0,npixs=0;
72   Float_t xyz[3];
73     
74   for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
75     {
76       Int_t prow = padrow - NRows[patch][0];
77       if(!track->GetCrossingPoint(padrow,xyz))  
78         {
79           printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
80           continue;
81         }
82       
83       fTransform->Slice2Sector(slice,padrow,sector,row);
84       fTransform->Local2Raw(xyz,sector,row);
85       npixs=0;
86       //Look at both sides of the crossing point:
87       for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
88         {
89           if(p<0 || p>fTransform->GetNPads(padrow)) continue;
90           lut_index = (prow<<8) + p;
91           if(track_lut[eta_index][lut_index]>0) //There was a signal here
92             npixs++;
93         }
94       
95       if(npixs > 0)
96         {
97           nrow++;
98         }         
99     }
100   if(nrow >= NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)//this was a good track
101     {
102       track->SetEtaIndex(eta_index);
103       if(remove)
104         RemoveTrackFromImage(track,eta_index);
105       return kTRUE;
106     }
107   else
108     return kFALSE;
109 }
110
111 void AliL3HoughEval::LookInsideRawRoad(AliL3TrackArray *tracks,Int_t eta_index,Bool_t remove)
112 {
113   //Evalaute the track candidates by looking along the trajectory.
114   //If remove is on, the pixels along the track will be removed. 
115
116
117   if(!fHoughTransformer)
118     {
119       printf("AliL3HoughEval: No transformer object\n");
120       return;
121     }
122   
123   Int_t patch = fHoughTransformer->fPatch;
124   Int_t slice = fHoughTransformer->fSlice;
125
126   
127   Int_t sector,row;
128   Int_t lut_index;
129   Char_t **track_lut = fHoughTransformer->fTrackTable;
130   Int_t nrow,npixs;
131   Float_t xyz[3];
132   
133   for(Int_t i=0; i<tracks->GetNTracks(); i++)
134     {
135       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
136       if(!track) {printf("No track\n"); return;}
137       
138       nrow = 0;
139       
140       for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
141         {
142           Int_t prow = padrow - NRows[patch][0];
143           if(!track->GetCrossingPoint(padrow,xyz))  
144             {
145               printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
146               continue;
147             }
148                   
149           fTransform->Slice2Sector(slice,padrow,sector,row);
150           fTransform->Local2Raw(xyz,sector,row);
151           npixs=0;
152           //Look at both sides of the crossing point:
153           for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
154             {
155               if(p<0 || p>fTransform->GetNPads(padrow)) continue;
156               lut_index = (prow<<8) + p;
157               if(track_lut[eta_index][lut_index]>0) //There was a signal here
158                 npixs++;
159             }
160           
161           if(npixs > 0)
162             {
163               nrow++;
164             }     
165         }
166       if(nrow < NRows[patch][1]-NRows[patch][0]-fNumOfRowsToMiss)
167         tracks->Remove(i); //this was not a good enough track
168       else if(remove)
169         RemoveTrackFromImage(track,eta_index); //this was a good track, so remove it from the image.
170     }
171   
172   tracks->Compress();
173 }
174
175 void AliL3HoughEval::RemoveTrackFromImage(AliL3HoughTrack *track,Int_t eta_index)
176 {
177   //Remove the pixels along the track in the image. 
178   
179   Int_t patch = fHoughTransformer->fPatch;
180   Int_t slice = fHoughTransformer->fSlice;
181
182   Int_t maxnpads = fTransform->GetNPads(NRows[patch][1]);
183   
184   Int_t lut_index;
185   Char_t **track_lut = fHoughTransformer->fTrackTable;
186   Int_t sector,row;
187   Float_t xyz[3];
188   for(Int_t padrow = NRows[patch][0]; padrow<=NRows[patch][1]; padrow++)
189     {
190       Int_t prow = padrow - NRows[patch][0];
191       if(!track->GetCrossingPoint(padrow,xyz))  
192         {
193           printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
194           continue;
195         }
196       
197       fTransform->Slice2Sector(slice,padrow,sector,row);
198       fTransform->Local2Raw(xyz,sector,row);
199       
200       for(Int_t p=(Int_t)xyz[1]-fNumOfPadsToLook; p<=(Int_t)xyz[1]+fNumOfPadsToLook; p++)
201         {
202           if(p<0 || p>fTransform->GetNPads(padrow)) continue;
203           lut_index = (prow<<8) + p;
204           track_lut[eta_index][lut_index] = -1; //remove it!!
205         }
206     }
207   
208   
209 }
210
211 void AliL3HoughEval::DisplaySlice(TH2F *hist)
212 {
213
214   AliL3Digits *pixel;
215   
216   for(Int_t padrow = 0; padrow < 174; padrow++)
217     {
218       
219       for(pixel=(AliL3Digits*)fHoughTransformer->fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
220         {
221           
222           Int_t sector,row;
223           Float_t xyz_pix[3];
224           fTransform->Slice2Sector(fHoughTransformer->fSlice,padrow,sector,row);
225           fTransform->Raw2Local(xyz_pix,sector,row,pixel->fPad,pixel->fTime); //y alone pad
226           
227           if(hist)
228             {
229               hist->Fill(xyz_pix[0],xyz_pix[1],pixel->fCharge);
230             }    
231           
232         }
233     }
234   
235 }
236
237 void AliL3HoughEval::DefineGoodParticles(Char_t *rootfile,Double_t pet)
238 {
239   //define the particles that produce good enough signals to be recognized in the transform
240
241   Int_t num_eta_segments = fHoughTransformer->fNumEtaSegments;
242   Double_t eta_max = fHoughTransformer->fEtaMax;
243   fMcTrackTable = new Int_t[num_eta_segments];
244   for(Int_t i=0; i<num_eta_segments; i++)
245     fMcTrackTable[i]=0;
246   Double_t etaslice = eta_max/num_eta_segments;
247
248   Int_t patch = fHoughTransformer->fPatch;
249   Int_t slice = fHoughTransformer->fSlice;
250   
251   TFile *file = new TFile(rootfile);
252   file->cd();
253
254   AliRun *gAlice = (AliRun*)file->Get("gAlice");
255   gAlice->GetEvent(0);
256   
257   TClonesArray *particles=gAlice->Particles();
258   Int_t np = particles->GetEntriesFast();
259   
260   Int_t row_to_miss = fNumOfRowsToMiss;
261   AliTPCParam *param = (AliTPCParam*)file->Get("75x40_100x60");
262   Int_t zero = param->GetZeroSup();
263   TTree *TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
264   AliSimDigits da, *digits=&da;
265   TD->GetBranch("Segment")->SetAddress(&digits); //Return pointer to branch segment.
266   Int_t *good=new Int_t[np];
267   Int_t *count = new Int_t[np]; //np number of particles.
268   Int_t good_number = NRows[patch][1]-NRows[patch][0]-row_to_miss;
269   Int_t i;
270   for (i=0; i<np; i++) count[i]=0;
271   Int_t sectors_by_rows=(Int_t)TD->GetEntries();
272   for (i=0; i<sectors_by_rows; i++) 
273     {
274       if (!TD->GetEvent(i)) continue;
275       Int_t sec,row;
276       param->AdjustSectorRow(digits->GetID(),sec,row);
277       Int_t ss,sr;
278       fTransform->Sector2Slice(ss,sr,sec,row);
279       if(ss!=slice) continue;
280       if(sr < NRows[patch][0]) continue;
281       if(sr > NRows[patch][1]) break;
282       digits->First();
283       while (digits->Next()) {
284         Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
285         Short_t dig = digits->GetDigit(it,ip);
286         Int_t idx0=digits->GetTrackID(it,ip,0); 
287         Int_t idx1=digits->GetTrackID(it,ip,1);
288         Int_t idx2=digits->GetTrackID(it,ip,2);
289         if (idx0>=0 && dig>=zero) count[idx0]+=1;
290         if (idx1>=0 && dig>=zero) count[idx1]+=1;
291         if (idx2>=0 && dig>=zero) count[idx2]+=1;
292       }
293       for (Int_t j=0; j<np; j++) 
294         {
295           if (count[j]>1) {//at least two digits at this padrow
296             good[j]++;
297           }
298           count[j]=0;
299         }
300     }
301   delete[] count;
302   
303   Int_t good_one=0;
304   //TObjArray *part = new TObjArray(0,0);
305   for(i=0; i<np; i++)
306    {
307      TParticle *p = (TParticle*)particles->UncheckedAt(i);
308      if(p->GetFirstMother()>0) continue; //secondary particle
309      if(good[i] < good_number) continue; //too few padrows
310      Double_t ptg=p->Pt();
311      if(ptg<pet) continue;
312      Int_t eta_index = (Int_t)(p->Eta()/etaslice);
313      if(eta_index < 0 || eta_index > num_eta_segments)
314        continue;
315      //fMcTrackTable[eta_index]++;
316      //part->AddLast(p);
317      good_one++;
318    }
319   
320   printf("nparticles %d\n",good_one);
321   file->Close();
322   delete [] good;
323   delete file;
324   //return part;
325 }
326
327
328 void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
329 {
330   
331   Int_t slice = fHoughTransformer->fSlice;
332   
333   TFile *file = new TFile(rootfile);
334   file->cd();
335   
336   AliRun *gAlice = (AliRun*)file->Get("gAlice");
337   gAlice->GetEvent(0);
338   
339   TClonesArray *particles=gAlice->Particles();  
340   Int_t n=particles->GetEntriesFast();
341   Float_t torad=TMath::Pi()/180;
342   Float_t phi_min = slice*20 - 10;
343   Float_t phi_max = slice*20 + 10;
344   
345   for (Int_t j=0; j<n; j++) {
346     TParticle *p=(TParticle*)particles->UncheckedAt(j);
347     if (p->GetFirstMother()>=0) continue;  //secondary particle
348     if(p->Eta() < eta[0] || p->Eta() > eta[1]) continue;
349     Double_t ptg=p->Pt(),pxg=p->Px(),pyg=p->Py();//,pzg=p->Pz();
350     Double_t phi_part = TMath::ATan2(pyg,pxg);
351     if (phi_part < 0) phi_part += 2*TMath::Pi();
352     
353     if(phi_part < phi_min*torad || phi_part > phi_max*torad) {continue;}
354     if(ptg<0.100) continue;
355     
356     AliL3HoughTrack *sel_track=0;
357     
358     Double_t min_dist = 10000;
359     
360     for(Int_t t=0; t<merged_tracks->GetNTracks(); t++)
361       {
362         AliL3HoughTrack *track = (AliL3HoughTrack*)merged_tracks->GetCheckedTrack(t);
363         if(!track) {printf("AliL3HoughEval: NO TRACK\n"); continue;}
364         Float_t phi0[2] = {track->GetPhi0(),0};
365         fTransform->Local2GlobalAngle(phi0,slice);
366         Double_t dPt = ptg - track->GetPt();
367         Double_t dPhi0 = phi_part - phi0[0];
368         Double_t distance = pow(dPt,2) + pow(dPhi0,2);
369         if(distance < min_dist)
370           {
371             min_dist = distance;
372             sel_track = track;
373           }      
374         
375       }
376     if(sel_track)
377       sel_track->SetBestMCid(j,min_dist);
378     
379     Double_t dpt = fabs(ptg-sel_track->GetPt());
380     Double_t dphi0 = fabs(phi_part-sel_track->GetPhi0());
381     //printf("Found match, min_dist %f dPt %f dPhi0 %f\n",min_dist,dpt,dphi0);
382   }
383   file->Close();
384   delete file;
385   
386 }