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 | |
21 | ClassImp(AliL3HoughEval) |
22 | |
23 | |
24 | AliL3HoughEval::AliL3HoughEval() |
25 | { |
26 | //Default constructor |
27 | fTransform = new AliL3Transform(); |
f80b98cb |
28 | fHoughTransformer = 0; |
29 | fNumOfRowsToMiss = 1; |
30 | fNumOfPadsToLook = 1; |
4de874d1 |
31 | } |
32 | |
33 | |
34 | AliL3HoughEval::AliL3HoughEval(AliL3HoughTransformer *transformer) |
35 | { |
36 | //Constructor |
37 | fHoughTransformer = transformer; |
38 | fTransform = new AliL3Transform(); |
f80b98cb |
39 | fNumOfRowsToMiss = 1; |
40 | fNumOfPadsToLook = 1; |
4de874d1 |
41 | } |
42 | |
43 | |
44 | AliL3HoughEval::~AliL3HoughEval() |
45 | { |
46 | //Destructor |
47 | if(fTransform) |
48 | delete fTransform; |
f80b98cb |
49 | if(fMcTrackTable) |
50 | delete [] fMcTrackTable; |
4de874d1 |
51 | } |
52 | |
53 | |
f80b98cb |
54 | Bool_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 | |
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 | |
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 | |
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; |
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 | |
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 | |
f80b98cb |
237 | void 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 | |
328 | void 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 | } |