Added several GET-Functions, fixed bug in old config.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformer.cxx
CommitLineData
b5a207b4 1//$Id$
2
b1886074 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ASV
4cafa5fc 5
48f3c46f 6#include "AliL3MemHandler.h"
e26acabd 7#include "AliL3Logging.h"
99e7186b 8#include "AliL3HoughTransformer.h"
4de874d1 9#include "AliL3Defs.h"
10#include "AliL3Transform.h"
4cafa5fc 11#include "AliL3DigitData.h"
99e7186b 12#include "AliL3Histogram.h"
4de874d1 13
b1886074 14//_____________________________________________________________
15// AliL3HoughTransformer
16//
17// Hough transformation class
18
4de874d1 19ClassImp(AliL3HoughTransformer)
20
21AliL3HoughTransformer::AliL3HoughTransformer()
22{
23 //Default constructor
4cafa5fc 24
4de874d1 25}
26
99e7186b 27AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments)
52a2a604 28{
52a2a604 29 fSlice = slice;
30 fPatch = patch;
99e7186b 31 fNEtaSegments = n_eta_segments;
32 fEtaMin = 0;
33 fEtaMax = fSlice < 18 ? 0.9 : -0.9;
34 fTransform = new AliL3Transform();
35 fThreshold = 0;
48f3c46f 36 fNDigitRowData=0;
37 fDigitRowData=0;
52a2a604 38}
39
4de874d1 40AliL3HoughTransformer::~AliL3HoughTransformer()
41{
4de874d1 42 if(fTransform)
43 delete fTransform;
99e7186b 44 DeleteHistograms();
4cafa5fc 45}
46
99e7186b 47void AliL3HoughTransformer::DeleteHistograms()
4cafa5fc 48{
99e7186b 49 if(!fParamSpace)
50 return;
51 for(Int_t i=0; i<fNEtaSegments; i++)
4fc9a6a4 52 {
53 if(!fParamSpace[i]) continue;
54 delete fParamSpace[i];
55 }
99e7186b 56 delete [] fParamSpace;
4cafa5fc 57}
58
e26acabd 59void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
60 Int_t nybin,Double_t phimin,Double_t phimax)
61{
62 //Set the minimum absolute pt value, and phi0 angles given in degrees.
63
e26acabd 64 Double_t bfact = 0.0029980;
65 Double_t bfield = 0.2;
66 Double_t x = bfact*bfield/pt_min;
b5a207b4 67 CreateHistograms(nxbin,-1.*x,x,nybin,phimin*ToRad,phimax*ToRad);
e26acabd 68}
69
99e7186b 70void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
71 Int_t nybin,Double_t ymin,Double_t ymax)
4cafa5fc 72{
4fc9a6a4 73
99e7186b 74 fParamSpace = new AliL3Histogram*[fNEtaSegments];
4cafa5fc 75
99e7186b 76 Char_t histname[256];
77 for(Int_t i=0; i<fNEtaSegments; i++)
4de874d1 78 {
99e7186b 79 sprintf(histname,"paramspace_%d",i);
80 fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
4de874d1 81 }
4cafa5fc 82}
83
e26acabd 84void AliL3HoughTransformer::Reset()
85{
48f3c46f 86 //Reset all the histograms
e26acabd 87
88 if(!fParamSpace)
89 {
90 LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
91 <<"No histograms to reset"<<ENDLOG;
92 return;
93 }
94
95 for(Int_t i=0; i<fNEtaSegments; i++)
96 fParamSpace[i]->Reset();
e26acabd 97}
98
99e7186b 99void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
9f33a1db 100{
99e7186b 101 fNDigitRowData = ndigits;
102 fDigitRowData = ptr;
4de874d1 103}
104
b5a207b4 105Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
106{
107
108 Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
109 Double_t index = (eta-fEtaMin)/etaslice;
110 return (Int_t)index;
111}
112
4fc9a6a4 113void AliL3HoughTransformer::TransformCircle()
4de874d1 114{
99e7186b 115 //Transform the input data with a circle HT.
116
52a2a604 117
99e7186b 118 //Set pointer to the data
119 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
120 if(!tempPt || fNDigitRowData==0)
4de874d1 121 {
4fc9a6a4 122 printf("\nAliL3HoughTransformer::TransformCircle : No input data!!!\n\n");
99e7186b 123 return;
4de874d1 124 }
99e7186b 125
b5a207b4 126 //Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
4cafa5fc 127 for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
128 {
99e7186b 129 AliL3DigitData *digPt = tempPt->fDigitData;
130 if(i != (Int_t)tempPt->fRow)
4cafa5fc 131 {
4fc9a6a4 132 printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
99e7186b 133 continue;
134 }
135 for(UInt_t j=0; j<tempPt->fNDigit; j++)
136 {
137 UShort_t charge = digPt[j].fCharge;
138 UChar_t pad = digPt[j].fPad;
139 UShort_t time = digPt[j].fTime;
48f3c46f 140 if(charge <= fThreshold)
99e7186b 141 continue;
4cafa5fc 142 Int_t sector,row;
99e7186b 143 Float_t xyz[3];
4cafa5fc 144 fTransform->Slice2Sector(fSlice,i,sector,row);
99e7186b 145 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
4fc9a6a4 146 Double_t eta = fTransform->GetEta(xyz);
b5a207b4 147 Int_t eta_index = GetEtaIndex(eta);//(Int_t)((eta-fEtaMin)/etaslice);
99e7186b 148 if(eta_index < 0 || eta_index >= fNEtaSegments)
149 continue;
4cafa5fc 150
e26acabd 151 //Get the correct histogrampointer:
99e7186b 152 AliL3Histogram *hist = fParamSpace[eta_index];
153 if(!hist)
4cafa5fc 154 {
4fc9a6a4 155 printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
99e7186b 156 continue;
4cafa5fc 157 }
4cafa5fc 158
99e7186b 159 //Start transformation
e26acabd 160 Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); // + xyz[2]*xyz[2]);
99e7186b 161 Float_t phi = fTransform->GetPhi(xyz);
4de874d1 162
99e7186b 163 //Fill the histogram along the phirange
164 for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
4de874d1 165 {
99e7186b 166 Float_t phi0 = hist->GetBinCenterY(b);
167 Float_t kappa = 2*sin(phi - phi0)/R;
168 hist->Fill(kappa,phi0,charge);
4de874d1 169 }
4de874d1 170 }
48f3c46f 171 AliL3MemHandler::UpdateRowPointer(tempPt);
4de874d1 172 }
4de874d1 173}
174
b1886074 175void AliL3HoughTransformer::TransformCircleC()
176{
177 //Circle transform, using combinations of every 2 points lying
178 //on different padrows and within the same etaslice.
179
180
181 Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
182 AliL3DigitRowData **rowPt = new AliL3DigitRowData*[nrows];
183
184 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
185 if(!tempPt)
186 printf("\nAliL3HoughTransformer::TransformCircleC() : Zero data pointer\n");
187
188 Int_t prow;
189 for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
190 {
191 prow = i - NRows[fPatch][0];
192 rowPt[prow] = tempPt;
193 AliL3MemHandler::UpdateRowPointer(tempPt);
194 }
b5a207b4 195 //Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
b1886074 196
197 AliL3DigitData *digPt;
198 Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
199 UShort_t charge1,charge2,time;
200 UChar_t pad;
201 Float_t xyz[3];
202 Int_t sector,row,eta_index1,eta_index2,tot_charge;
203 for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
204 {
205 prow = i - NRows[fPatch][0];
206 digPt = rowPt[prow]->fDigitData;
207 for(UInt_t di=0; di<rowPt[prow]->fNDigit; di++)
208 {
209 charge1 = digPt[di].fCharge;
210 pad = digPt[di].fPad;
211 time = digPt[di].fTime;
212 if(charge1 <= fThreshold)
213 continue;
214 fTransform->Slice2Sector(fSlice,i,sector,row);
215 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
216 eta = fTransform->GetEta(xyz);
b5a207b4 217 eta_index1 = GetEtaIndex(eta);//(Int_t)((eta-fEtaMin)/etaslice);
b1886074 218 if(eta_index1 < 0 || eta_index1 >= fNEtaSegments)
219 continue;
220 r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
221 phi1 = atan2(xyz[1],xyz[0]);
222
223 //Get the correct histogrampointer:
224 AliL3Histogram *hist = fParamSpace[eta_index1];
225 if(!hist)
226 {
227 printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",eta_index1);
228 continue;
229 }
230
231 for(Int_t j=i+1; j<NRows[fPatch][1]; j++)
232 {
233 prow = j - NRows[fPatch][0];
234 digPt = rowPt[prow]->fDigitData;
235 for(UInt_t ni=0; ni<rowPt[prow]->fNDigit; ni++)
236 {
237 charge2 = digPt[ni].fCharge;
238 pad = digPt[ni].fPad;
239 time = digPt[ni].fTime;
240 if(charge2 <= fThreshold)
241 continue;
242 fTransform->Slice2Sector(fSlice,j,sector,row);
243 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
244 eta = fTransform->GetEta(xyz);
b5a207b4 245 eta_index2 = GetEtaIndex(eta);//(Int_t)((eta-fEtaMin)/etaslice);
b1886074 246 if(eta_index2 != eta_index1)
247 continue;
248 r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
249 phi2 = atan2(xyz[1],xyz[0]);
250
251 phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
252 kappa = 2*sin(phi2-phi_0)/r2;
253 tot_charge = charge1+charge2;
254 //printf("Filling kappa %f phi %f charge %d\n",kappa,phi_0,tot_charge);
255 hist->Fill(kappa,phi_0,tot_charge);
256 }
257 }
258 }
259 }
260
261 delete [] rowPt;
262}
263
4fc9a6a4 264void AliL3HoughTransformer::TransformLine()
265{
266 //Do a line transform on the data.
267
268
269 AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
270 if(!tempPt || fNDigitRowData==0)
271 {
272 printf("\nAliL3HoughTransformer::TransformLine : No input data!!!\n\n");
273 return;
274 }
275
b5a207b4 276 //Double_t etaslice = (fEtaMax - fEtaMin)/fNEtaSegments;
4fc9a6a4 277 for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
278 {
279 AliL3DigitData *digPt = tempPt->fDigitData;
280 if(i != (Int_t)tempPt->fRow)
281 {
282 printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
283 continue;
284 }
285 for(UInt_t j=0; j<tempPt->fNDigit; j++)
286 {
287 UShort_t charge = digPt[j].fCharge;
288 UChar_t pad = digPt[j].fPad;
289 UShort_t time = digPt[j].fTime;
290 if(charge < fThreshold)
291 continue;
292 Int_t sector,row;
293 Float_t xyz[3];
294 fTransform->Slice2Sector(fSlice,i,sector,row);
295 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
296 Float_t eta = fTransform->GetEta(xyz);
b5a207b4 297 Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
4fc9a6a4 298 if(eta_index < 0 || eta_index >= fNEtaSegments)
299 continue;
300
301 //Get the correct histogram:
302 AliL3Histogram *hist = fParamSpace[eta_index];
303 if(!hist)
304 {
305 printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
306 continue;
307 }
308 for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
309 {
310 Double_t theta = hist->GetBinCenterX(xbin);
311 Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
312 hist->Fill(theta,rho,charge);
313 }
314 }
48f3c46f 315 AliL3MemHandler::UpdateRowPointer(tempPt);
4fc9a6a4 316 }
317
318}