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