This commit was generated by cvs2svn to compensate for changes in r3174,
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformer.cxx
CommitLineData
4de874d1 1#include <math.h>
2
3#include <TFile.h>
4#include <TTree.h>
5#include <TH2.h>
6
7#include "AliTPCParam.h"
8#include "AliSimDigits.h"
9
10#include "AliL3Defs.h"
11#include "AliL3Transform.h"
12#include "AliL3HoughPixel.h"
13#include "AliL3HoughTransformer.h"
14#include "AliL3HoughTrack.h"
15#include "AliL3TrackArray.h"
16
17ClassImp(AliL3HoughTransformer)
18
19AliL3HoughTransformer::AliL3HoughTransformer()
20{
21 //Default constructor
22 fTransform = 0;
23 fEtaMin = 0;
24 fEtaMax = 0;
25 fSlice = 0;
26 fPatch = 0;
27 fRowContainer = 0;
28 fPhiRowContainer = 0;
29 fNumOfPadRows=0;
30 fContainerBounds=0;
31 fNDigits=0;
32 fIndex = 0;
33}
34
35
36AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments)
37{
38 //Constructor
39
40 fTransform = new AliL3Transform();
41
42
43 fEtaMin = etarange[0];
44 fEtaMax = etarange[1];
45 fSlice = slice;
46 fPatch = patch;
47 fNPhiSegments = phi_segments;
48 fNumOfPadRows=NRowsSlice;
49}
50
51
52AliL3HoughTransformer::~AliL3HoughTransformer()
53{
54 //Destructor
55 if(fRowContainer)
56 delete [] fRowContainer;
57 if(fTransform)
58 delete fTransform;
59 if(fPhiRowContainer)
60 delete [] fPhiRowContainer;
61 if(fIndex)
62 {
63 for(Int_t i=0; i<fNDigits; i++)
64 delete [] fIndex[i];
65 delete [] fIndex;
66 }
67}
68
69void AliL3HoughTransformer::InitTemplates(TH2F *hist)
70{
71
72 AliL3Digits *pixel;
73
74 Int_t ymin = hist->GetYaxis()->GetFirst();
75 Int_t ymax = hist->GetYaxis()->GetLast();
76 Int_t nbinsy = hist->GetNbinsY();
77
78 fIndex = new Int_t*[fNDigits];
79 for(Int_t i=0; i<fNDigits; i++)
80 fIndex[i] = new Int_t[nbinsy+1];
81
82 Int_t sector,row;
83 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
84 {
85
86 for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
87 {
88 Float_t xyz[3];
89 fTransform->Slice2Sector(fSlice,padrow,sector,row);
90 fTransform->Raw2Local(xyz,sector,row,pixel->fPad,pixel->fTime);
91
92 Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
93
94 Double_t phi_pix = fTransform->GetPhi(xyz);
95 Short_t signal = pixel->fCharge;
96 Int_t index = pixel->fIndex;
97 if(index >= fNDigits)
98 printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
99 for(Int_t p=ymin; p<=ymax; p++)
100 {
101
102 Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
103 Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
104
105 Int_t bin = hist->FindBin(kappa,phi0);
106 if(fIndex[index][p]!=0)
107 printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
108 fIndex[index][p] = bin;
109 }
110 }
111 }
112
113}
114
115
116void AliL3HoughTransformer::CountBins()
117{
118
119 Int_t middle_row = 87; //middle of the slice
120
121 Double_t r_in_bundle = fTransform->Row2X(middle_row);
122 // Double_t phi_min = (fSlice*20 - 10)*ToRad;
123 //Double_t phi_max = (fSlice*20 + 10)*ToRad;
124 Double_t phi_min = -15*ToRad;
125 Double_t phi_max = 15*ToRad;
126
127 Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
128 Double_t min_phi0 = 10000;
129 Double_t max_phi0 = 0;
130 Double_t min_kappa = 100000;
131 Double_t max_kappa = 0;
132
133 Int_t xbin = 60;
134 Int_t ybin = 60;
135 Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
136 Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
137
138 TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
139
140 for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
141 {
142 for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
143 {
144 for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
145 {
146 Float_t xyz[3];
147 Int_t sector,row;
148 fTransform->Slice2Sector(fSlice,padrow,sector,row);
149 fTransform->Raw2Global(xyz,sector,row,pad,time);
150 Double_t eta = fTransform->GetEta(xyz);
151 if(eta < fEtaMin || eta > fEtaMax) continue;
152 fTransform->Global2Local(xyz,sector);
153 Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
154 Double_t phi_pix = fTransform->GetPhi(xyz);
155
156 for(Int_t p=0; p<=fNPhiSegments; p++)
157 {
158 Double_t phi_in_bundle = phi_min + p*phi_slice;
159
160 Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
161
162 Double_t phi0 = atan(tanPhi0);
163 // if(phi0 < 0.55 || phi0 > 0.88) continue;
164
165 //if(phi0 < 0) phi0 = phi0 +2*Pi;
166 //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
167
168 Double_t angle = phi_pix - phi0;
169 Double_t kappa = 2*sin(angle)/r_pix;
170 histo->Fill(kappa,phi0,1);
171 if(phi0 < min_phi0)
172 min_phi0 = phi0;
173 if(phi0 > max_phi0)
174 max_phi0 = phi0;
175 if(kappa < min_kappa)
176 min_kappa = kappa;
177 if(kappa > max_kappa)
178 max_kappa = kappa;
179
180 }
181
182 }
183
184 }
185
186 }
187 Int_t count=0,bi=0;
188
189 Int_t xmin = histo->GetXaxis()->GetFirst();
190 Int_t xmax = histo->GetXaxis()->GetLast();
191 Int_t ymin = histo->GetYaxis()->GetFirst();
192 Int_t ymax = histo->GetYaxis()->GetLast();
193
194
195 for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
196 {
197 for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
198 {
199 bi++;
200 Int_t bin = histo->GetBin(xbin,ybin);
201 if(histo->GetBinContent(bin)>0)
202 count++;
203 }
204 }
205
206
207 printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
208 printf("Phi, min %f max %f\n",min_phi0,max_phi0);
209 printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
210 histo->Draw("box");
211}
212
213
214void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row)
215{
216 //Transformation is done with respect to local coordinates in slice.
217
218
219 AliL3Digits *pix1;
220 Int_t sector,row;
221
222 //Define a common point
223
224 /*
225 Double_t rowdist1 = fTransform->Row2X(middle_row-1);
226 Double_t rowdist2 = fTransform->Row2X(middle_row);
227 Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2;
228 */
229
230 Double_t r_in_bundle = fTransform->Row2X(middle_row);
231 //Make overlap between slices
232 Double_t phi_min = -15*ToRad;
233 Double_t phi_max = 15*ToRad;
234
235 Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
236
237 for(Int_t p=0; p <= fNPhiSegments; p++)
238 {
239 Double_t phi_in_bundle = phi_min + p*phi_slice;
240 //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle);
241
242 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
243 //for(Int_t padrow = middle_row; padrow <= 173; padrow++)
244 //for(Int_t padrow = 0; padrow <= middle_row; padrow++)
245 {
246
247 for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
248 {
249
250 Float_t xyz[3];
251 fTransform->Slice2Sector(fSlice,padrow,sector,row);
252 fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
253 //fTransform->Raw2Global(xyz,sector,row,pix1->fPad,pix1->fTime);
254
255 Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
256
257 Double_t phi_pix = fTransform->GetPhi(xyz);
258
259 //Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
260 Double_t tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle));
261
262 Double_t phi0 = atan(tanPhi0);
263 //if(padrow > middle_row)
264 //phi0 = -phi0;
265 //if(phi0 < 0.55 || phi0 > 0.88) continue;
266
267 Double_t angle = phi_pix - phi0;
268 Double_t kappa = 2*sin(angle)/r_pix;
269
270 //Double_t angle = phi_in_bundle - phi0;
271 //Double_t kappa = 2*sin(angle)/r_in_bundle;
272
273 //if(kappa < -0.006 || kappa > 0.006) continue;
274
275 Short_t signal = pix1->fCharge;
276
277 hist->Fill(kappa,phi0,signal);
278
279
280 }
281
282 }
283 }
284
285}
286
287
288void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
289{
290 //Transformation is done with respect to local coordinates in slice.
291 //Transform every pixel into whole phirange, using parametrisation:
292 //kappa = 2*sin(phi-phi0)/R
293
294 printf("Transforming 1 pixel only\n");
295
296 AliL3Digits *pix1;
297 Int_t sector,row;
298
299 Int_t ymin = hist->GetYaxis()->GetFirst();
300 Int_t ymax = hist->GetYaxis()->GetLast();
301 Int_t nbinsy = hist->GetNbinsY();
302
303 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
304 {
305
306 for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
307 {
308
309 Short_t signal = pix1->fCharge;
310 Int_t index = pix1->fIndex;
311
312 for(Int_t p=0; p <= nbinsy; p++)
313 hist->AddBinContent(fIndex[index][p],signal);
314
315 }
316
317 }
318
319
320}
321
322/*
323 void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
324 {
325 //Transformation is done with respect to local coordinates in slice.
326 //Transform every pixel into whole phirange, using parametrisation:
327 //kappa = 2*sin(phi-phi0)/R
328
329 printf("Transforming 1 pixel only\n");
330
331 AliL3Digits *pix1;
332 Int_t sector,row;
333
334 Int_t ymin = hist->GetYaxis()->GetFirst();
335 Int_t ymax = hist->GetYaxis()->GetLast();
336
337 for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
338 {
339
340 for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
341 {
342
343 Float_t xyz[3];
344 fTransform->Slice2Sector(fSlice,padrow,sector,row);
345 fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
346
347 Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
348
349 Double_t phi_pix = fTransform->GetPhi(xyz);
350 Short_t signal = pix1->fCharge;
351
352 for(Int_t p=ymin+1; p<=ymax; p++)
353 {
354 Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
355 Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
356 hist->Fill(kappa,phi0,signal);
357 }
358
359 }
360
361 }
362
363
364 }
365*/
366
367void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
368{
369
370 for(Int_t i=0; i<tracks->GetNTracks(); i++)
371 {
372 AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
373 if(!track) {printf("AliL3HoughTransformer::TransformLines2Circle : NO TRACK IN ARRAY\n"); continue;}
374
375 Double_t xmin = fTransform->Row2X(track->GetFirstRow());
376 Double_t xmax = fTransform->Row2X(track->GetLastRow());
377
378 Double_t a = -1./tan(track->GetPsiLine());
379 Double_t b = track->GetDLine()/sin(track->GetPsiLine());
380
381 Double_t ymin = a*xmin + b;
382 Double_t ymax = a*xmax + b;
383
384 Double_t middle_x = xmin + (xmax-xmin)/2;
385 Double_t middle_y = ymin + (ymax-ymin)/2;
386
387 Double_t r_middle = sqrt(middle_x*middle_x + middle_y*middle_y);
388 Double_t phi = atan2(middle_y,middle_x);
389 Double_t phi0 = 2*phi - track->GetPsiLine();
390 Double_t kappa = 2*sin(phi-phi0)/r_middle;
391 hist->Fill(kappa,phi0,track->GetWeight());
392
393 }
394
395
396}
397
398void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw)
399{
400 //Transform every pixel into histogram, using parametrisation:
401 //D = x*cos(psi) + y*sin(psi)
402
403 // printf("In Transform; rowrange %d %d ref_row %d phirange %f %f\n",rowrange[0],rowrange[1],ref_row,phirange[0],phirange[1]);
404
405 AliL3Digits *pix1;
406
407 Int_t xmin = hist->GetXaxis()->GetFirst();
408 Int_t xmax = hist->GetXaxis()->GetLast();
409
410 Double_t x0 = fTransform->Row2X(ref_row);
411 Double_t y0 = 0;
412
413 Int_t sector,row;
414 //for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
415
416 Double_t phi_min = -10*ToRad;
417 Double_t phi_max = 10*ToRad;
418 Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
419
420
421 Int_t phi_min_index = (Int_t)((phirange[0]*ToRad-phi_min)/delta_phi);
422 Int_t phi_max_index = (Int_t)((phirange[1]*ToRad-phi_min)/delta_phi);
423
424
425 Int_t index;
426
427 for(Int_t phi=phi_min_index; phi <= phi_max_index; phi++)
428 {
429 for(Int_t padrow = rowrange[0]; padrow <= rowrange[1]; padrow++)
430 {
431
432 index = phi*fNumOfPadRows + padrow;
433 //printf("Looping index %d\n",index);
434 if(index > fContainerBounds || index < 0)
435 {
436 printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
437 return;
438 }
439 //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
440 for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
441 {
442 //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
443 Float_t xyz[3];
444 fTransform->Slice2Sector(fSlice,padrow,sector,row);
445 fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
446
447 if(raw)
448 raw->Fill(xyz[0],xyz[1],pix1->fCharge);
449
450 xyz[0] = xyz[0]-x0;
451 xyz[1] = xyz[1]-y0;
452
453 //printf("calculating...");
454 for(Int_t d=xmin+1; d<=xmax; d++)
455 {
456 Double_t psi = hist->GetXaxis()->GetBinCenter(d);
457 Double_t D = xyz[0]*cos(psi) + xyz[1]*sin(psi);
458
459 Short_t signal = pix1->fCharge;
460 hist->Fill(psi,D,signal);
461 //printf("Filling histo, psi %f D %f\n",psi,D);
462 }
463 //printf("done\n");
464 }
465 //printf(" done\n");
466 }
467
468 }
469
470}
471
472void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
473{
474
475 TFile *file = new TFile(rootfile);
476 file->cd();
477
478 AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60");
479 TTree *t=(TTree*)file->Get("TreeD_75x40_100x60");
480
481 AliSimDigits da, *digarr=&da;
482 t->GetBranch("Segment")->SetAddress(&digarr);
483 Stat_t num_of_entries=t->GetEntries();
484
485 Int_t digit_counter=0;
486 Float_t xyz[3];
487 Double_t eta,phi;
488
489 Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
490 printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
491
492 if(fRowContainer)
493 delete [] fRowContainer;
494 fRowContainer = new AliL3HoughContainer[fNumOfPadRows];
495 memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
496
497
498 fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
499 printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
500 if(fPhiRowContainer)
501 delete [] fPhiRowContainer;
502 fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
503 memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
504
505 Double_t phi_min = -10*ToRad;
506 Double_t phi_max = 10*ToRad;
507 Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
508 Int_t index;
509 digit_counter=0;
510
511 for (Int_t i=0; i<num_of_entries; i++)
512 {
513 t->GetEvent(i);
514 Int_t sector;
515 Int_t row;
516 param->AdjustSectorRow(digarr->GetID(),sector,row);
517 Int_t slice,padrow;
518 fTransform->Sector2Slice(slice,padrow,sector,row);
519 if(slice != fSlice) continue;
520 if(padrow < NRows[fPatch][0]) continue;
521 if(padrow > NRows[fPatch][1]) break;
522 digarr->First();
523 do {
524 Int_t time=digarr->CurrentRow();
525 Int_t pad=digarr->CurrentColumn();
526 Short_t signal=digarr->CurrentDigit();
527 if(time < param->GetMaxTBin()-1 && time > 0)
528 if(digarr->GetDigit(time-1,pad) <= param->GetZeroSup()
529 && digarr->GetDigit(time+1,pad) <= param->GetZeroSup())
530 continue;
531
532
533 fTransform->Raw2Global(xyz,sector,row,pad,time);
534 eta = fTransform->GetEta(xyz);
535 if(eta < fEtaMin || eta > fEtaMax) continue;
536 fTransform->Global2Local(xyz,sector);
537
538
539 phi = fTransform->GetPhi(xyz);
540 if(hist)
541 hist->Fill(xyz[0],xyz[1],signal);
542
543 AliL3Digits *dig = new AliL3Digits;
544 dig->fIndex = digit_counter;
545 digit_counter++;
546 dig->fCharge = signal;
547 dig->fPad = pad;
548 dig->fTime = time;
549
550 if(fRowContainer[padrow].first == NULL)
551 fRowContainer[padrow].first = (void*)dig;
552 else
553 ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
554 fRowContainer[padrow].last = (void*)dig;
555
556 Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
557 index = phi_index*fNumOfPadRows + padrow;
558 if(phi_index > fContainerBounds || phi_index < 0)
559 {
560 printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
561 continue;
562 }
563
564 if(fPhiRowContainer[index].first == NULL)
565 fPhiRowContainer[index].first = (void*)dig;
566 else
567 ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
568 fPhiRowContainer[index].last=(void*)dig;
569
570 }while (digarr->Next());
571
572 }
573
574 fNDigits = digit_counter;
575 printf("digitcounter %d\n",digit_counter);
576 printf("Allocated %d bytes to pixels\n",digit_counter*sizeof(AliL3Digits));
577 file->Close();
578 delete file;
579
580}