]>
Commit | Line | Data |
---|---|---|
d5d754e5 | 1 | //$Id$ |
2 | ||
3 | // Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de> | |
4 | //*-- Copyright & copy CL | |
5 | ||
6 | #include "AliL3StandardIncludes.h" | |
7 | ||
8 | #include "AliL3Logging.h" | |
9 | #include "AliL3HoughTransformerLUT.h" | |
10 | #include "AliL3Transform.h" | |
11 | #include "AliL3MemHandler.h" | |
12 | #include "AliL3DigitData.h" | |
13 | #include "AliL3Histogram.h" | |
14 | ||
15 | #if GCCVERSION == 3 | |
16 | using namespace std; | |
17 | #endif | |
18 | ||
19 | //_____________________________________________________________ | |
20 | // AliL3HoughTransformerLUT | |
21 | // | |
22 | // Hough transformation class using Look-UP-Tables | |
23 | // | |
24 | ||
25 | ClassImp(AliL3HoughTransformerLUT) | |
26 | ||
27 | AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer() | |
28 | { | |
29 | fMinRow=0; | |
30 | fMaxRow=0; | |
31 | ||
32 | fNRows=0; | |
33 | fNEtas=0; | |
34 | fNPhi0=0; | |
35 | fSector=0; | |
36 | fSectorRow=0; | |
37 | fZSign=0; | |
38 | fZLengthPlusOff=0; | |
39 | fTimeWidth=0; | |
40 | fPadPitch=0; | |
41 | fEtaSlice=0; | |
42 | ||
43 | fLUTX=0; | |
44 | fLUTY=0; | |
45 | fLUTEta=0; | |
46 | fLUTphi0=0; | |
47 | fLUT2sinphi0=0; | |
48 | fLUT2cosphi0=0; | |
49 | } | |
50 | ||
51 | AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments) | |
52 | { | |
53 | AliL3HoughTransformerLUT(); | |
54 | ||
55 | Init(slice,patch,n_eta_segments); | |
56 | } | |
57 | ||
58 | AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT() | |
59 | { | |
60 | DeleteHistograms(); | |
61 | ||
62 | #ifdef do_mc | |
63 | if(fTrackID) | |
64 | { | |
65 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
66 | { | |
67 | if(!fTrackID[i]) continue; | |
68 | delete fTrackID[i]; | |
69 | } | |
70 | delete [] fTrackID; | |
71 | } | |
72 | #endif | |
73 | ||
74 | if(fNRows){ | |
75 | delete[] fLUTX; | |
76 | delete[] fLUTY; | |
77 | fNRows=0; | |
78 | } | |
79 | if(fNEtas){ | |
80 | delete[] fLUTEta; | |
81 | fNEtas=0; | |
82 | } | |
83 | } | |
84 | ||
85 | void AliL3HoughTransformerLUT::DeleteHistograms() | |
86 | { | |
87 | if(fNPhi0){ | |
88 | delete[] fLUT2sinphi0; | |
89 | delete[] fLUT2cosphi0; | |
90 | delete[] fLUTphi0; | |
91 | fNPhi0=0; | |
92 | fLUTphi0=0; | |
93 | fLUT2sinphi0=0; | |
94 | fLUT2cosphi0=0; | |
95 | } | |
96 | ||
97 | if(!fParamSpace){ | |
98 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
99 | { | |
100 | if(!fParamSpace[i]) continue; | |
101 | delete fParamSpace[i]; | |
102 | } | |
103 | delete [] fParamSpace; | |
104 | } | |
105 | } | |
106 | ||
107 | void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments) | |
108 | { | |
109 | AliL3HoughBaseTransformer::Init(slice,patch,n_eta_segments); | |
110 | ||
111 | //delete old LUT tables | |
112 | if(fNRows){ | |
113 | fNRows=0; | |
114 | delete[] fLUTX; | |
115 | delete[] fLUTY; | |
116 | } | |
117 | if(fNEtas){ | |
118 | delete[] fLUTEta; | |
119 | fNEtas=0; | |
120 | } | |
121 | ||
122 | //member values | |
123 | fMinRow=AliL3Transform::GetFirstRow(patch);; | |
124 | fMaxRow=AliL3Transform::GetLastRow(patch);; | |
125 | fNRows=AliL3Transform::GetNRows(patch);; | |
126 | fNEtas=n_eta_segments; | |
127 | ||
128 | AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow); | |
129 | fZSign = slice < 18 ? 1:-1; | |
130 | fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset(); | |
131 | fTimeWidth=AliL3Transform::GetZWidth(); | |
132 | ||
133 | fPadPitch=0; | |
134 | if(fSector<AliL3Transform::GetNSectorLow()) | |
135 | fPadPitch=AliL3Transform::GetPadPitchWidthLow(); | |
136 | else | |
137 | fPadPitch=AliL3Transform::GetPadPitchWidthUp(); | |
138 | ||
139 | Float_t etamax_=GetEtaMax(); | |
140 | Float_t etamin_=GetEtaMin(); | |
141 | fEtaSlice=(etamax_-etamin_)/n_eta_segments; | |
142 | ||
143 | //lookup tables for X and Y | |
144 | fLUTX=new Float_t[fNRows]; | |
145 | fLUTY=new Float_t[fNRows]; | |
146 | for(Int_t rr=0;rr<fNRows;rr++){ | |
147 | fLUTX[rr]=Float_t(AliL3Transform::Row2X(rr+fMinRow)); | |
148 | fLUTY[rr]=Float_t(0.5*(AliL3Transform::GetNPads(rr+fMinRow)-1)*fPadPitch); | |
149 | //cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl; | |
150 | } | |
151 | ||
152 | //lookup tables for rz2s <=> etas | |
153 | fLUTEta=new Float_t[fNEtas]; | |
154 | for(Int_t rr=0;rr<fNEtas;rr++){ | |
155 | fLUTEta[rr]=CalcRoverZ2(etamin_+(rr+1)*fEtaSlice); | |
156 | //cout << rr << ": " << fLUTEta[rr] << endl; | |
157 | } | |
158 | } | |
159 | ||
160 | void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax) | |
161 | { | |
162 | //Create the histograms (parameter space). | |
163 | //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis). | |
164 | //The arguments give the range and binning; | |
165 | //nxbin = #bins in kappa | |
166 | //nybin = #bins in phi0 | |
167 | //pt_min = mimium Pt of track (corresponding to maximum kappa) | |
168 | //phi_min = mimimum phi0 (degrees) | |
169 | //phi_max = maximum phi0 (degrees) | |
170 | ||
171 | Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min; | |
172 | Double_t torad = AliL3Transform::Pi()/180; | |
173 | CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad); | |
174 | } | |
175 | ||
176 | void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax) | |
177 | { | |
178 | fParamSpace = new AliL3Histogram*[fNEtas]; | |
179 | ||
180 | Char_t histname[256]; | |
181 | for(Int_t i=0; i<fNEtas; i++) | |
182 | { | |
183 | sprintf(histname,"paramspace_%d",i); | |
184 | fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax); | |
185 | } | |
186 | ||
187 | #ifdef do_mc | |
188 | if(fDoMC) | |
189 | { | |
190 | AliL3Histogram *hist = fParamSpace[0]; | |
191 | Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2); | |
192 | cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl; | |
193 | fTrackID = new TrackIndex*[GetNEtaSegments()]; | |
194 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
195 | fTrackID[i] = new TrackIndex[ncells]; | |
196 | } | |
197 | #endif | |
198 | ||
199 | //create lookup table for sin and cos | |
200 | fNPhi0=nybin+1; | |
201 | ||
202 | fLUTphi0=new Float_t[fNPhi0]; | |
203 | fLUT2sinphi0=new Float_t[fNPhi0]; | |
204 | fLUT2cosphi0=new Float_t[fNPhi0]; | |
205 | Float_t diff=(ymax-ymin)/nybin; | |
206 | Float_t phi0=ymin-0.5*diff; | |
207 | for(Int_t i=0; i<fNPhi0; i++){ | |
208 | phi0+=diff; | |
209 | fLUTphi0[i]=phi0; | |
210 | fLUT2sinphi0[i]=Float_t(2*sin(phi0)); | |
211 | fLUT2cosphi0[i]=Float_t(2*cos(phi0)); | |
212 | cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl; | |
213 | } | |
214 | } | |
215 | ||
216 | void AliL3HoughTransformerLUT::Reset() | |
217 | { | |
218 | //Reset all the histograms | |
219 | ||
220 | if(!fParamSpace) | |
221 | { | |
222 | LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms") | |
223 | <<"No histograms to reset"<<ENDLOG; | |
224 | return; | |
225 | } | |
226 | ||
227 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
228 | fParamSpace[i]->Reset(); | |
229 | #ifdef do_mc | |
230 | if(fDoMC) | |
231 | { | |
232 | AliL3Histogram *hist = fParamSpace[0]; | |
233 | Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2); | |
234 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
235 | memset(fTrackID[i],0,ncells*sizeof(TrackIndex)); | |
236 | } | |
237 | #endif | |
238 | } | |
239 | ||
240 | ||
241 | Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta) | |
242 | { | |
243 | //Return the histogram index of the corresponding eta. | |
244 | Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments(); | |
245 | Double_t index = (eta-GetEtaMin())/etaslice; | |
246 | return (Int_t)index; | |
247 | } | |
248 | ||
249 | inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index) | |
250 | { | |
251 | if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0) | |
252 | return 0; | |
253 | if(!fParamSpace[eta_index]) | |
254 | return 0; | |
255 | return fParamSpace[eta_index]; | |
256 | } | |
257 | ||
258 | Float_t AliL3HoughTransformerLUT::CalcRoverZ2(Float_t eta) | |
259 | { | |
260 | Float_t e=exp(2*eta); | |
261 | Float_t ret=(e+1)/(e-1); | |
262 | ret*=ret; | |
263 | return ret; | |
264 | } | |
265 | ||
266 | Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2) | |
267 | { | |
268 | Float_t rz=sqrt(roverz2); | |
269 | if(fZSign>0) rz=-rz; | |
270 | Float_t ret=(1+rz)/(rz-1); | |
271 | ret=0.5*log(ret); | |
272 | return ret; | |
273 | } | |
274 | ||
275 | Float_t AliL3HoughTransformerLUT::CalcX(Int_t row) | |
276 | { | |
277 | return fLUTX[row]; | |
278 | } | |
279 | ||
280 | Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row) | |
281 | { | |
282 | return pad*fPadPitch-fLUTY[row]; | |
283 | } | |
284 | ||
285 | Float_t AliL3HoughTransformerLUT::CalcZ(Int_t time) | |
286 | { | |
287 | Float_t ret=time*fTimeWidth; | |
288 | if(fZSign>0) ret=fZLengthPlusOff-ret; | |
289 | else ret=ret-fZLengthPlusOff; | |
290 | return ret; | |
291 | } | |
292 | ||
293 | Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice) | |
294 | { | |
295 | Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments(); | |
296 | Double_t eta=(Double_t)((eta_index+0.5)*eta_slice); | |
297 | if(slice>17) eta*=-1; | |
298 | return eta; | |
299 | } | |
300 | ||
301 | void AliL3HoughTransformerLUT::TransformCircle() | |
302 | { | |
303 | #if 0 | |
304 | //Transform the input data with a circle HT. | |
305 | //The function loops over all the data, and transforms each pixel with the equations: | |
306 | // | |
307 | //kappa = 2/R*sin(phi - phi0) | |
308 | // | |
309 | //where R = sqrt(x*x +y*y), and phi = arctan(y/x) | |
310 | // | |
311 | //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find | |
312 | //which histogram in which the pixel should be transformed, the eta-value is calcluated | |
313 | //and the proper histogram index is found by GetEtaIndex(eta). | |
314 | ||
315 | ||
316 | AliL3DigitRowData *tempPt = GetDataPointer(); | |
317 | if(!tempPt) | |
318 | { | |
319 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data") | |
320 | <<"No input data "<<ENDLOG; | |
321 | return; | |
322 | } | |
323 | ||
324 | //Loop over the padrows: | |
325 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) | |
326 | { | |
327 | //Get the data on this padrow: | |
328 | AliL3DigitData *digPt = tempPt->fDigitData; | |
329 | if(i != (Int_t)tempPt->fRow) | |
330 | { | |
331 | cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl; | |
332 | continue; | |
333 | } | |
334 | ||
335 | //Loop over the data on this padrow: | |
336 | for(UInt_t j=0; j<tempPt->fNDigit; j++) | |
337 | { | |
338 | UShort_t charge = digPt[j].fCharge; | |
339 | UChar_t pad = digPt[j].fPad; | |
340 | UShort_t time = digPt[j].fTime; | |
341 | if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold()) | |
342 | continue; | |
343 | Int_t sector,row; | |
344 | Float_t xyz[3]; | |
345 | ||
346 | //Transform data to local cartesian coordinates: | |
347 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); | |
348 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
349 | ||
350 | //Calculate the eta: | |
351 | Double_t eta = AliL3Transform::GetEta(xyz); | |
352 | ||
353 | //Get the corresponding index, which determines which histogram to fill: | |
354 | Int_t eta_index = GetEtaIndex(eta); | |
355 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) | |
356 | continue; | |
357 | ||
358 | //Get the correct histogrampointer: | |
359 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
360 | if(!hist) | |
361 | { | |
362 | printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index); | |
363 | continue; | |
364 | } | |
365 | ||
366 | //Do the transformation: | |
367 | Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
368 | Float_t phi = AliL3Transform::GetPhi(xyz); | |
369 | ||
370 | //Fill the histogram along the phirange | |
371 | for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++) | |
372 | { | |
373 | Float_t phi0 = hist->GetBinCenterY(b); | |
374 | Float_t kappa = 2*sin(phi - phi0)/R; | |
375 | hist->Fill(kappa,phi0,charge); | |
376 | #ifdef do_mc | |
377 | if(fDoMC) | |
378 | { | |
379 | Int_t bin = hist->FindBin(kappa,phi0); | |
380 | for(Int_t t=0; t<3; t++) | |
381 | { | |
382 | Int_t label = digPt[j].fTrackID[t]; | |
383 | if(label < 0) break; | |
384 | UInt_t c; | |
385 | for(c=0; c<MaxTrack; c++) | |
386 | if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0) | |
387 | break; | |
388 | if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl; | |
389 | fTrackID[eta_index][bin].fLabel[c] = label; | |
390 | fTrackID[eta_index][bin].fNHits[c]++; | |
391 | } | |
392 | } | |
393 | #endif | |
394 | } | |
395 | } | |
396 | ||
397 | //Move the data pointer to the next padrow: | |
398 | AliL3MemHandler::UpdateRowPointer(tempPt); | |
399 | } | |
400 | #endif //to do later | |
401 | } | |
402 | ||
403 | Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi) | |
404 | { | |
405 | if(!fDoMC) | |
406 | { | |
407 | cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl; | |
408 | return -1; | |
409 | } | |
410 | ||
411 | #ifdef do_mc | |
412 | if(eta_index < 0 || eta_index > GetNEtaSegments()) | |
413 | { | |
414 | cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl; | |
415 | return -1; | |
416 | } | |
417 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
418 | Int_t bin = hist->FindBin(kappa,psi); | |
419 | Int_t label=-1; | |
420 | Int_t max=0; | |
421 | for(UInt_t i=0; i<MaxTrack; i++) | |
422 | { | |
423 | Int_t nhits=fTrackID[eta_index][bin].fNHits[i]; | |
424 | if(nhits == 0) break; | |
425 | if(nhits > max) | |
426 | { | |
427 | max = nhits; | |
428 | label = fTrackID[eta_index][bin].fLabel[i]; | |
429 | } | |
430 | } | |
431 | return label; | |
432 | #endif | |
433 | cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl; | |
434 | return -1; | |
435 | } | |
436 | ||
437 | void AliL3HoughTransformerLUT::Print() | |
438 | { | |
439 | cout << "fSlice: " << GetSlice() << endl; | |
440 | cout << "fPatch: " << GetPatch() << endl; | |
441 | cout << "fSector: " << fSector << endl; | |
442 | cout << "fSectorRow: " << fSectorRow << endl; | |
443 | cout << "fMinRow: " << fMinRow << endl; | |
444 | cout << "fMaxRow: " << fMaxRow << endl; | |
445 | cout << "fNRows: " << fNRows << endl; | |
446 | cout << "fNEtas: " << fNEtas << endl; | |
447 | cout << "fNPhi0: " << fNPhi0 << endl; | |
448 | cout << "fZSign: " << fZSign << endl; | |
449 | cout << "fZLengthPlusOff: " << fZLengthPlusOff << endl; | |
450 | cout << "fPadPitch: " << fPadPitch << endl; | |
451 | cout << "fTimeWidth: " << fTimeWidth << endl; | |
452 | ||
453 | if(!fNRows) return; | |
454 | cout << "fLUTX " << fNRows << endl; | |
455 | for(Int_t i=0;i<fNRows;i++) cout << "fLUTX[" << i << "]=" << (Float_t)fLUTX[i] << endl; | |
456 | cout << "fLUTY " << fNRows << endl; | |
457 | for(Int_t i=0;i<fNRows;i++) cout << "fLUTY[" << i << "]=" << fLUTY[i] << endl; | |
458 | if(!fNEtas) return; | |
459 | cout << "fLUTEta " << fNEtas << endl; | |
460 | for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl; | |
461 | if(!fNPhi0) return; | |
462 | cout << "fLUTphi0 " << fNPhi0 << endl; | |
463 | for(Int_t i=0;i<fNPhi0;i++) cout << "fLUTPhi0[" << i << "]=" << fLUTphi0[i] << endl; | |
464 | cout << "fLUT2sinphi0 " << fNPhi0 << endl; | |
465 | for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2sinphi0[" << i << "]=" << fLUT2sinphi0[i] << endl; | |
466 | cout << "fLUT2cosphi0 " << fNPhi0 << endl; | |
467 | for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2cosphi0[" << i << "]=" << fLUT2cosphi0[i] << endl; | |
468 | } |