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