La Navale continue ses frappes
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughEval.cxx
CommitLineData
f80b98cb 1//Author: Anders Strand Vestbo
2//Last Modified: 28.6.01
3
4de874d1 4#include <TClonesArray.h>
5#include <TH2.h>
6#include <TH1.h>
7#include <TFile.h>
8#include <AliRun.h>
9#include <TParticle.h>
f80b98cb 10#include <TTree.h>
4de874d1 11
f80b98cb 12#include "AliTPCParam.h"
13#include "AliSimDigits.h"
4de874d1 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
21ClassImp(AliL3HoughEval)
22
23
24AliL3HoughEval::AliL3HoughEval()
25{
26 //Default constructor
27 fTransform = new AliL3Transform();
f80b98cb 28 fHoughTransformer = 0;
29 fNumOfRowsToMiss = 1;
30 fNumOfPadsToLook = 1;
4de874d1 31}
32
33
34AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer)
35{
36 //Constructor
37 fHoughTransformer = transformer;
38 fTransform = new AliL3Transform();
f80b98cb 39 fNumOfRowsToMiss = 1;
40 fNumOfPadsToLook = 1;
4de874d1 41}
42
43
44AliL3HoughEval::~AliL3HoughEval()
45{
46 //Destructor
47 if(fTransform)
48 delete fTransform;
f80b98cb 49 if(fMcTrackTable)
50 delete [] fMcTrackTable;
4de874d1 51}
52
53
f80b98cb 54Bool_t AliL3HoughEval::LookInsideRoad(AliL3HoughTrack *track,Int_t eta_index,Bool_t remove)
4de874d1 55{
56 //Look at rawdata along the road specified by the track candidates.
f80b98cb 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
111void 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
4de874d1 116
117 if(!fHoughTransformer)
118 {
119 printf("AliL3HoughEval: No transformer object\n");
120 return;
121 }
122
4de874d1 123 Int_t patch = fHoughTransformer->fPatch;
124 Int_t slice = fHoughTransformer->fSlice;
4de874d1 125
f80b98cb 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
4de874d1 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
f80b98cb 138 nrow = 0;
4de874d1 139
140 for(Int_t padrow = NRows[patch][0]; padrow <= NRows[patch][1]; padrow++)
141 {
f80b98cb 142 Int_t prow = padrow - NRows[patch][0];
4de874d1 143 if(!track->GetCrossingPoint(padrow,xyz))
144 {
145 printf("AliL3HoughEval::LookInsideRoad : Track does not cross line!!\n");
146 continue;
147 }
f80b98cb 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++)
4de874d1 154 {
f80b98cb 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++;
4de874d1 159 }
f80b98cb 160
4de874d1 161 if(npixs > 0)
162 {
163 nrow++;
164 }
165 }
f80b98cb 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.
4de874d1 170 }
171
172 tracks->Compress();
f80b98cb 173}
174
175void 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;
4de874d1 181
f80b98cb 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
4de874d1 208
209}
210
211void 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
f80b98cb 237void AliL3HoughEval::DefineGoodParticles(Char_t *rootfile,Double_t pet)
4de874d1 238{
f80b98cb 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;
4de874d1 250
4de874d1 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();
f80b98cb 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);
9ee44abe 321 file->Close();
f80b98cb 322 delete [] good;
9ee44abe 323 delete file;
f80b98cb 324 //return part;
9ee44abe 325}
326
327
328void AliL3HoughEval::CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta)
329{
330
f80b98cb 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);
9ee44abe 338
f80b98cb 339 TClonesArray *particles=gAlice->Particles();
4de874d1 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 }
f80b98cb 383 file->Close();
384 delete file;
385
4de874d1 386}