Rewriting and cleaning up
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformer.cxx
1 //Author:        Anders Strand Vestbo
2 //Last Modified: 14.09.01
3
4 #include "AliL3HoughTransformer.h"
5 #include "AliL3Defs.h"
6 #include "AliL3Transform.h"
7 #include "AliL3DigitData.h"
8 #include "AliL3Histogram.h"
9
10 ClassImp(AliL3HoughTransformer)
11
12 AliL3HoughTransformer::AliL3HoughTransformer()
13 {
14   //Default constructor
15   
16 }
17
18 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments)
19 {
20   fSlice = slice;
21   fPatch = patch;
22   fNEtaSegments = n_eta_segments;
23   fEtaMin = 0;
24   fEtaMax = fSlice < 18 ? 0.9 : -0.9;
25   fTransform = new AliL3Transform();
26   fThreshold = 0;
27   
28 }
29
30 AliL3HoughTransformer::~AliL3HoughTransformer()
31 {
32   if(fTransform)
33     delete fTransform;
34   DeleteHistograms();
35 }
36
37 void AliL3HoughTransformer::DeleteHistograms()
38 {
39   if(!fParamSpace)
40     return;
41   for(Int_t i=0; i<fNEtaSegments; i++)
42     {
43       if(!fParamSpace[i]) continue;
44       delete fParamSpace[i];
45     }
46   delete [] fParamSpace;
47 }
48
49 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
50                                              Int_t nybin,Double_t ymin,Double_t ymax)
51 {
52   
53   fParamSpace = new AliL3Histogram*[fNEtaSegments];
54   
55   Char_t histname[256];
56   for(Int_t i=0; i<fNEtaSegments; i++)
57     {
58       sprintf(histname,"paramspace_%d",i);
59       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
60     }
61 }
62
63 void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
64 {
65   fNDigitRowData = ndigits;
66   fDigitRowData = ptr;
67 }
68
69 void AliL3HoughTransformer::UpdateDataPointer(AliL3DigitRowData *&tempPt)
70 {
71   //Update the data pointer to the next padrow
72   
73   Byte_t *tmp = (Byte_t*)tempPt;
74   Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
75   tmp += size;
76   tempPt = (AliL3DigitRowData*)tmp;
77 }
78
79 void AliL3HoughTransformer::TransformCircle()
80 {
81   //Transform the input data with a circle HT.
82   
83
84   //Set pointer to the data
85   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
86   if(!tempPt || fNDigitRowData==0)
87     {
88       printf("\nAliL3HoughTransformer::TransformCircle : No input data!!!\n\n");
89       return;
90     }
91   
92   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
93   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
94     {
95       AliL3DigitData *digPt = tempPt->fDigitData;
96       if(i != (Int_t)tempPt->fRow)
97         {
98           printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
99           continue;
100         }
101       for(UInt_t j=0; j<tempPt->fNDigit; j++)
102         {
103           UShort_t charge = digPt[j].fCharge;
104           UChar_t pad = digPt[j].fPad;
105           UShort_t time = digPt[j].fTime;
106           if(charge < fThreshold)
107             continue;
108           Int_t sector,row;
109           Float_t xyz[3];
110           fTransform->Slice2Sector(fSlice,i,sector,row);
111           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
112           Double_t eta = fTransform->GetEta(xyz);
113           Int_t eta_index = (Int_t)(eta/etaslice);
114           if(eta_index < 0 || eta_index >= fNEtaSegments)
115             continue;
116           
117           //Get the correct histogram:
118           AliL3Histogram *hist = fParamSpace[eta_index];
119           if(!hist)
120             {
121               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
122               continue;
123             }
124
125           //Start transformation
126           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1] + xyz[2]*xyz[2]);
127           Float_t phi = fTransform->GetPhi(xyz);
128           
129           //Fill the histogram along the phirange
130           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
131             {
132               Float_t phi0 = hist->GetBinCenterY(b);
133               Float_t kappa = 2*sin(phi - phi0)/R;
134               hist->Fill(kappa,phi0,charge);
135             }
136         }
137       UpdateDataPointer(tempPt);
138     }
139 }
140
141 void AliL3HoughTransformer::TransformLine()
142 {
143   //Do a line transform on the data.
144
145   
146   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
147   if(!tempPt || fNDigitRowData==0)
148     {
149       printf("\nAliL3HoughTransformer::TransformLine : No input data!!!\n\n");
150       return;
151     }
152   
153   Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
154   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
155     {
156       AliL3DigitData *digPt = tempPt->fDigitData;
157       if(i != (Int_t)tempPt->fRow)
158         {
159           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
160           continue;
161         }
162       for(UInt_t j=0; j<tempPt->fNDigit; j++)
163         {
164           UShort_t charge = digPt[j].fCharge;
165           UChar_t pad = digPt[j].fPad;
166           UShort_t time = digPt[j].fTime;
167           if(charge < fThreshold)
168             continue;
169           Int_t sector,row;
170           Float_t xyz[3];
171           fTransform->Slice2Sector(fSlice,i,sector,row);
172           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
173           Float_t eta = fTransform->GetEta(xyz);
174           Int_t eta_index = (Int_t)(eta/etaslice);
175           if(eta_index < 0 || eta_index >= fNEtaSegments)
176             continue;
177           
178           //Get the correct histogram:
179           AliL3Histogram *hist = fParamSpace[eta_index];
180           if(!hist)
181             {
182               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
183               continue;
184             }
185           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
186             {
187               Double_t theta = hist->GetBinCenterX(xbin);
188               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
189               hist->Fill(theta,rho,charge);
190             }
191         }
192       UpdateDataPointer(tempPt);
193     }
194   
195 }