Removed obsolete code, and removed a typo
[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 "AliL3HoughTransformer.h"
13 #include "AliL3HoughTrack.h"
14 #include "AliL3TrackArray.h"
15
16 ClassImp(AliL3HoughTransformer)
17
18 AliL3HoughTransformer::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;
28   fVolume=0;
29   fNumOfPadRows=0;
30   fContainerBounds=0;
31   fNDigits=0;
32   fIndex = 0;
33 }
34
35
36 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange)
37 {
38   //Constructor
39   
40   fTransform = new AliL3Transform();
41
42   fEtaMin = etarange[0];
43   fEtaMax = etarange[1];
44   fSlice = slice; 
45   fPatch = patch;
46   fNumOfPadRows=NRowsSlice;
47 }
48
49 AliL3HoughTransformer::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
71
72 AliL3HoughTransformer::~AliL3HoughTransformer()
73 {
74   //Destructor
75   if(fRowContainer)
76     delete [] fRowContainer;
77   if(fTransform)
78     delete fTransform;
79   if(fPhiRowContainer)
80     delete [] fPhiRowContainer;
81   if(fVolume)
82     delete [] fVolume;
83   if(fIndex)
84     {
85       for(Int_t i=0; i<fNDigits; i++)
86         delete [] fIndex[i];
87       delete [] fIndex;
88     }
89 }
90
91 void 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;
105   
106   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
107     {        
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);
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;
126               //printf("kappa %f phi0 %f\n",kappa,phi0);
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
138 void 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
236 void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
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
241   //Assumes you run InitTemplates firstly!!!!
242
243   AliL3Digits *pix1;
244     
245   Int_t nbinsy = hist->GetNbinsY();
246
247   Int_t totsignal=0,vol_index;
248   if(fNumEtaSegments==1)
249     eta_index=0; //only looking in one etaslice.
250
251   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
252     {
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)
256         {
257           
258           Short_t signal = pix1->fCharge;
259           Int_t index = pix1->fIndex;
260           if(index < 0) continue; //This pixel has been removed.
261           totsignal += signal;
262           for(Int_t p=0; p <= nbinsy; p++)
263             hist->AddBinContent(fIndex[index][p],signal);
264                           
265         }
266       
267     }
268     
269   printf("Total signal %d\n",totsignal);
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
317 void 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
348 void 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             }
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)
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
422 void 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];
437   Double_t eta;
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
448   //fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
449   //printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
450
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;
468   Int_t index;
469   digit_counter=0;
470
471   printf("\nLoading ALL pixels in slice\n\n");
472
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);
497         
498         if(eta < fEtaMin || eta > fEtaMax) continue;
499         fTransform->Global2Local(xyz,sector);
500         
501         
502         //phi = fTransform->GetPhi(xyz);
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         
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)
523           {
524             
525             //printf("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
526             continue;
527           }
528         
529         if(fVolume[index].first == NULL)
530           fVolume[index].first = (void*)dig;
531         else
532           ((AliL3Digits*)(fVolume[index].last))->fNextVolumePixel = dig;
533         fVolume[index].last = (void*)dig;
534         
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         */
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 }