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