]>
Commit | Line | Data |
---|---|---|
28ac1a53 | 1 | //$Id$ |
2 | ||
3 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> | |
4 | //*-- Copyright © ASV | |
5 | ||
e06900d5 | 6 | #include "AliL3StandardIncludes.h" |
7 | ||
28ac1a53 | 8 | #include "AliL3Logging.h" |
9 | #include "AliL3HoughClusterTransformer.h" | |
e06900d5 | 10 | #include "AliL3MemHandler.h" |
28ac1a53 | 11 | #include "AliL3SpacePointData.h" |
12 | #include "AliL3Transform.h" | |
13 | #include "AliL3DigitData.h" | |
14 | #include "AliL3Histogram.h" | |
15 | #include "AliL3ClustFinderNew.h" | |
16 | ||
e06900d5 | 17 | #if GCCVERSION == 3 |
18 | using namespace std; | |
19 | #endif | |
20 | ||
28ac1a53 | 21 | //_____________________________________________________________ |
22 | // AliL3HoughClusterTransformer | |
23 | // | |
24 | // Hough transformation class. | |
25 | // | |
26 | // This class performs the Hough transform on _clusters_, and not raw data. | |
27 | // Therefore, it first finds the clusters using the HLT cluster finder, and | |
28 | // then uses these as input for the transform. | |
29 | ||
30 | ClassImp(AliL3HoughClusterTransformer) | |
31 | ||
32 | AliL3HoughClusterTransformer::AliL3HoughClusterTransformer() | |
33 | { | |
34 | //Default constructor | |
35 | fParamSpace = 0; | |
36 | fNClusters = 0; | |
37 | fClusters = 0; | |
38 | fMemHandler = 0; | |
39 | #ifdef do_mc | |
40 | fTrackID = 0; | |
41 | #endif | |
42 | } | |
43 | ||
44 | AliL3HoughClusterTransformer::AliL3HoughClusterTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments) | |
45 | { | |
46 | //Normal constructor | |
47 | fParamSpace = 0; | |
48 | fNClusters=0; | |
49 | fClusters=0; | |
50 | fMemHandler=0; | |
51 | #ifdef do_mc | |
52 | fTrackID=0; | |
53 | #endif | |
54 | } | |
55 | ||
56 | AliL3HoughClusterTransformer::~AliL3HoughClusterTransformer() | |
57 | { | |
58 | DeleteHistograms(); | |
59 | if(fMemHandler) | |
60 | delete fMemHandler; | |
61 | #ifdef do_mc | |
62 | if(fTrackID) | |
63 | { | |
64 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
65 | { | |
66 | if(!fTrackID[i]) continue; | |
67 | delete fTrackID[i]; | |
68 | } | |
69 | delete [] fTrackID; | |
70 | } | |
71 | #endif | |
72 | } | |
73 | ||
74 | //void AliL3HoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){} | |
75 | ||
76 | void AliL3HoughClusterTransformer::DeleteHistograms() | |
77 | { | |
78 | if(fParamSpace) | |
79 | { | |
80 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
81 | { | |
82 | if(!fParamSpace[i]) continue; | |
83 | delete fParamSpace[i]; | |
84 | } | |
85 | delete [] fParamSpace; | |
86 | } | |
87 | } | |
88 | ||
89 | void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min, | |
90 | Int_t nybin,Double_t phimin,Double_t phimax) | |
91 | { | |
92 | //Create the histograms (parameter space). | |
93 | //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis). | |
94 | //The arguments give the range and binning; | |
95 | //nxbin = #bins in kappa | |
96 | //nybin = #bins in phi0 | |
97 | //pt_min = mimium Pt of track (corresponding to maximum kappa) | |
98 | //phi_min = mimimum phi0 (degrees) | |
99 | //phi_max = maximum phi0 (degrees) | |
100 | ||
101 | Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min; | |
102 | Double_t torad = AliL3Transform::Pi()/180; | |
103 | CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad); | |
104 | } | |
105 | ||
106 | void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax, | |
107 | Int_t nybin,Double_t ymin,Double_t ymax) | |
108 | { | |
109 | Char_t histname[256]; | |
110 | ||
111 | fParamSpace = new AliL3Histogram*[GetNEtaSegments()]; | |
112 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
113 | { | |
114 | sprintf(histname,"paramspace_%d",i); | |
115 | fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax); | |
116 | } | |
117 | #ifdef do_mc | |
118 | Int_t ncells = (nxbin+2)*(nybin+2); | |
119 | cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl; | |
120 | fTrackID = new TrackIndex*[GetNEtaSegments()]; | |
121 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
122 | fTrackID[i] = new TrackIndex[ncells]; | |
123 | #endif | |
124 | ||
125 | } | |
126 | ||
127 | void AliL3HoughClusterTransformer::Reset() | |
128 | { | |
129 | //Reset all the histograms | |
130 | ||
131 | if(!fParamSpace) | |
132 | { | |
133 | LOG(AliL3Log::kWarning,"AliL3HoughClusterTransformer::Reset","Histograms") | |
134 | <<"No histograms to reset"<<ENDLOG; | |
135 | return; | |
136 | } | |
137 | ||
138 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
139 | fParamSpace[i]->Reset(); | |
140 | ||
141 | if(fMemHandler) | |
142 | delete fMemHandler; | |
143 | fNClusters=0; | |
144 | fClusters=0; | |
145 | #ifdef do_mc | |
146 | AliL3Histogram *hist = fParamSpace[0]; | |
147 | Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2); | |
148 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
149 | memset(fTrackID[i],0,ncells*sizeof(TrackIndex)); | |
150 | #endif | |
151 | } | |
152 | ||
153 | Int_t AliL3HoughClusterTransformer::GetEtaIndex(Double_t eta) | |
154 | { | |
155 | //Return the histogram index of the corresponding eta. | |
156 | ||
157 | Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments(); | |
158 | Double_t index = (eta-GetEtaMin())/etaslice; | |
159 | return (Int_t)index; | |
160 | } | |
161 | ||
162 | inline AliL3Histogram *AliL3HoughClusterTransformer::GetHistogram(Int_t eta_index) | |
163 | { | |
164 | if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0) | |
165 | return 0; | |
166 | if(!fParamSpace[eta_index]) | |
167 | return 0; | |
168 | return fParamSpace[eta_index]; | |
169 | } | |
170 | ||
171 | Double_t AliL3HoughClusterTransformer::GetEta(Int_t eta_index,Int_t slice) | |
172 | { | |
173 | Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments(); | |
174 | Double_t eta=(Double_t)((eta_index+0.5)*eta_slice); | |
175 | if(slice>17) eta*=-1; | |
176 | return eta; | |
177 | } | |
178 | ||
179 | void AliL3HoughClusterTransformer::FindClusters() | |
180 | { | |
181 | //Find the clusters | |
182 | if(!GetDataPointer()) | |
183 | { | |
184 | cerr<<"AliL3HoughClusterTransformer::FindClusters : Zero data pointer"<<endl; | |
185 | return; | |
186 | } | |
187 | ||
188 | const Int_t maxpoints=100000; | |
189 | const Int_t pointsize = maxpoints * sizeof(AliL3SpacePointData); | |
190 | ||
191 | fMemHandler = new AliL3MemHandler(); | |
192 | fClusters = (AliL3SpacePointData*)fMemHandler->Allocate(pointsize); | |
193 | AliL3ClustFinderNew *cf = new AliL3ClustFinderNew(); | |
194 | cf->InitSlice(0,GetPatch(),AliL3Transform::GetFirstRow(GetPatch()),AliL3Transform::GetLastRow(GetPatch()),maxpoints); | |
195 | cf->SetDeconv(kFALSE); | |
196 | cf->SetOutputArray(fClusters); | |
197 | cf->Read(1,GetDataPointer()); | |
198 | cf->ProcessDigits(); | |
199 | fNClusters = cf->GetNumberOfClusters(); | |
200 | delete cf; | |
201 | } | |
202 | ||
203 | void AliL3HoughClusterTransformer::TransformCircle() | |
204 | { | |
205 | //Transform the input data with a circle HT. | |
206 | //The function loops over all the data, and transforms each cluster with the equations: | |
207 | // | |
208 | //kappa = 2/R*sin(phi - phi0) | |
209 | // | |
210 | //where R = sqrt(x*x +y*y), and phi = arctan(y/x) | |
211 | // | |
212 | //Each cluster then transforms into a curve in the (kappa,phi0)-space. In order to find | |
213 | //which histogram in which the pixel should be transformed, the eta-value is calcluated | |
214 | //and the proper histogram index is found by GetEtaIndex(eta). | |
215 | ||
216 | FindClusters(); | |
217 | if(!fClusters) | |
218 | { | |
219 | LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircle","Data") | |
220 | <<"No input data "<<ENDLOG; | |
221 | return; | |
222 | } | |
223 | ||
224 | Float_t xyz[3]; | |
225 | for(Int_t i=0; i<fNClusters; i++) | |
226 | { | |
227 | xyz[0] = fClusters[i].fX; | |
228 | xyz[1] = fClusters[i].fY; | |
229 | xyz[2] = fClusters[i].fZ; | |
230 | Double_t eta = AliL3Transform::GetEta(xyz); | |
231 | Int_t eta_index = GetEtaIndex(eta); | |
232 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) | |
233 | continue; | |
234 | ||
235 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
236 | ||
237 | //Do the transformation: | |
238 | Double_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
239 | Double_t phi = AliL3Transform::GetPhi(xyz); | |
240 | ||
241 | //Fill the histogram along the phirange | |
242 | for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++) | |
243 | { | |
244 | Double_t phi0 = hist->GetBinCenterY(b); | |
245 | Double_t kappa = 2*sin(phi - phi0)/R; | |
246 | hist->Fill(kappa,phi0,1); | |
247 | #ifdef do_mc | |
248 | Int_t bin = hist->FindBin(kappa,phi0); | |
249 | for(Int_t t=0; t<3; t++) | |
250 | { | |
251 | Int_t label = fClusters[i].fTrackID[t]; | |
252 | if(label < 0) break; | |
253 | UInt_t c; | |
254 | for(c=0; c<MaxTrack; c++) | |
255 | if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0) | |
256 | break; | |
257 | if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl; | |
258 | fTrackID[eta_index][bin].fLabel[c] = label; | |
259 | fTrackID[eta_index][bin].fNHits[c]++; | |
260 | } | |
261 | #endif | |
262 | } | |
263 | } | |
264 | } | |
265 | ||
266 | void AliL3HoughClusterTransformer::TransformCircleC(Int_t row_range) | |
267 | { | |
268 | //Circle transform, using combinations of every 2 points lying | |
269 | //on different padrows and within the same etaslice. | |
270 | ||
271 | FindClusters(); | |
272 | if(!fClusters) | |
273 | LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircleC","Data") | |
274 | <<"No input data "<<ENDLOG; | |
275 | ||
276 | Float_t xyz[3]; | |
277 | Double_t eta,r1,phi1,r2,phi2,kappa,psi; | |
278 | Int_t index1,index2; | |
279 | for(Int_t i=0; i<fNClusters; i++) | |
280 | { | |
281 | xyz[0] = fClusters[i].fX; | |
282 | xyz[1] = fClusters[i].fY; | |
283 | xyz[2] = fClusters[i].fZ; | |
284 | r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
285 | phi1 = atan2(xyz[1],xyz[0]); | |
286 | eta = AliL3Transform::GetEta(xyz); | |
287 | index1 = GetEtaIndex(eta); | |
288 | if(index1 < 0 || index1 >= GetNEtaSegments()) continue; | |
289 | ||
290 | for(Int_t j=i+1; j<fNClusters; j++) | |
291 | { | |
292 | if(fClusters[j].fPadRow == fClusters[i].fPadRow) continue; | |
293 | xyz[0] = fClusters[j].fX; | |
294 | xyz[1] = fClusters[j].fY; | |
295 | xyz[2] = fClusters[j].fZ; | |
296 | eta = AliL3Transform::GetEta(xyz); | |
297 | index2 = GetEtaIndex(eta); | |
298 | if(index1 != index2) continue; | |
299 | ||
300 | AliL3Histogram *hist = fParamSpace[index1]; | |
301 | ||
302 | r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
303 | phi2 = atan2(xyz[1],xyz[0]); | |
304 | ||
305 | psi = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) ); | |
306 | kappa = 2*sin(phi1-psi) / r1; | |
307 | hist->Fill(kappa,psi,1); | |
308 | #ifdef do_mc | |
309 | Int_t bin = hist->FindBin(kappa,psi); | |
310 | for(Int_t l=0; l<3; l++) | |
311 | { | |
312 | for(Int_t m=0; m<3; m++) | |
313 | { | |
314 | if(fClusters[i].fTrackID[l] == fClusters[j].fTrackID[m]) | |
315 | { | |
316 | Int_t label = fClusters[i].fTrackID[l]; | |
317 | if(label < 0) continue; | |
318 | UInt_t c; | |
319 | for(c=0; c<MaxTrack; c++) | |
320 | if(fTrackID[index1][bin].fLabel[c] == label || fTrackID[index1][bin].fNHits[c] == 0) | |
321 | break; | |
322 | if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl; | |
323 | fTrackID[index1][bin].fLabel[c] = label; | |
324 | fTrackID[index1][bin].fNHits[c]++; | |
325 | } | |
326 | } | |
327 | } | |
328 | #endif | |
329 | } | |
330 | } | |
331 | } | |
332 | ||
333 | ||
334 | Int_t AliL3HoughClusterTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi) | |
335 | { | |
336 | ||
337 | #ifdef do_mc | |
338 | if(eta_index < 0 || eta_index > GetNEtaSegments()) | |
339 | { | |
340 | cerr<<"AliL3HoughClusterTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl; | |
341 | return -1; | |
342 | } | |
343 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
344 | Int_t bin = hist->FindBin(kappa,psi); | |
345 | Int_t label=-1; | |
346 | Int_t max=0; | |
347 | for(UInt_t i=0; i<MaxTrack; i++) | |
348 | { | |
349 | Int_t nhits=fTrackID[eta_index][bin].fNHits[i]; | |
350 | if(nhits == 0) break; | |
351 | if(nhits > max) | |
352 | { | |
353 | max = nhits; | |
354 | label = fTrackID[eta_index][bin].fLabel[i]; | |
355 | } | |
356 | } | |
357 | return label; | |
358 | #endif | |
359 | cout<<"AliL3HoughClusterTransformer::GetTrackID : Compile with do_mc flag!"<<endl; | |
360 | return -1; | |
361 | } | |
362 |