]>
Commit | Line | Data |
---|---|---|
3e87ef69 | 1 | // @(#) $Id$ |
b5a207b4 | 2 | |
b1886074 | 3 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> |
3e87ef69 | 4 | //*-- Copyright © ALICE HLT Group |
4cafa5fc | 5 | |
e06900d5 | 6 | #include "AliL3StandardIncludes.h" |
7 | ||
e26acabd | 8 | #include "AliL3Logging.h" |
99e7186b | 9 | #include "AliL3HoughTransformer.h" |
e06900d5 | 10 | #include "AliL3MemHandler.h" |
4de874d1 | 11 | #include "AliL3Transform.h" |
4cafa5fc | 12 | #include "AliL3DigitData.h" |
93657878 | 13 | #include "AliL3HistogramAdaptive.h" |
4de874d1 | 14 | |
0bd0c1ef | 15 | #if __GNUC__ == 3 |
e06900d5 | 16 | using namespace std; |
17 | #endif | |
18 | ||
1f1942b8 | 19 | /** \class AliL3HoughTransformer |
20 | <pre> | |
b1886074 | 21 | //_____________________________________________________________ |
22 | // AliL3HoughTransformer | |
23 | // | |
24 | // Hough transformation class | |
3ef466c5 | 25 | // |
1f1942b8 | 26 | </pre> |
27 | */ | |
b1886074 | 28 | |
4de874d1 | 29 | ClassImp(AliL3HoughTransformer) |
30 | ||
31 | AliL3HoughTransformer::AliL3HoughTransformer() | |
32 | { | |
33 | //Default constructor | |
c52cf5d8 | 34 | fParamSpace = 0; |
44c7f8de | 35 | fDoMC = kFALSE;; |
3e87ef69 | 36 | fEtaOverlap=kFALSE; |
1eef4efc | 37 | #ifdef do_mc |
38 | fTrackID = 0; | |
39 | #endif | |
4de874d1 | 40 | } |
41 | ||
3e87ef69 | 42 | AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoEtaOverlap,Bool_t DoMC) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments) |
52a2a604 | 43 | { |
3ef466c5 | 44 | //Normal constructor |
237d3f5c | 45 | fParamSpace = 0; |
3e87ef69 | 46 | fDoMC = kFALSE; |
47 | fEtaOverlap = DoEtaOverlap; | |
48 | fDoMC=kFALSE; | |
1eef4efc | 49 | #ifdef do_mc |
50 | fTrackID = 0; | |
51 | #endif | |
52a2a604 | 52 | } |
53 | ||
4de874d1 | 54 | AliL3HoughTransformer::~AliL3HoughTransformer() |
55 | { | |
99e7186b | 56 | DeleteHistograms(); |
1eef4efc | 57 | #ifdef do_mc |
58 | if(fTrackID) | |
59 | { | |
60 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
61 | { | |
62 | if(!fTrackID[i]) continue; | |
63 | delete fTrackID[i]; | |
64 | } | |
65 | delete [] fTrackID; | |
66 | } | |
67 | #endif | |
4cafa5fc | 68 | } |
69 | ||
99e7186b | 70 | void AliL3HoughTransformer::DeleteHistograms() |
4cafa5fc | 71 | { |
99e7186b | 72 | if(!fParamSpace) |
73 | return; | |
c52cf5d8 | 74 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
4fc9a6a4 | 75 | { |
76 | if(!fParamSpace[i]) continue; | |
77 | delete fParamSpace[i]; | |
78 | } | |
99e7186b | 79 | delete [] fParamSpace; |
4cafa5fc | 80 | } |
81 | ||
b2a02bce | 82 | void AliL3HoughTransformer::CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres, |
83 | Int_t nybin,Float_t psi) | |
84 | { | |
85 | //Create histograms. | |
86 | //_Only_ to be used in case of the adaptive histograms! | |
87 | //phimax is given in radians!! | |
88 | ||
89 | if(ptmin > ptmax) | |
90 | { | |
1f1942b8 | 91 | cerr<<"AliL3HoughTransformer::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl; |
b2a02bce | 92 | return; |
93 | } | |
94 | if(psi < 0) | |
95 | { | |
96 | cerr<<"AliL3HoughTransformer::CreateHistograms: Wrong psi-angle "<<psi<<endl; | |
97 | return; | |
98 | } | |
99 | ||
100 | fParamSpace = new AliL3Histogram*[GetNEtaSegments()]; | |
101 | Char_t histname[256]; | |
102 | Int_t i; | |
103 | for(i=0; i<GetNEtaSegments(); i++) | |
104 | { | |
105 | sprintf(histname,"paramspace_%d",i); | |
106 | fParamSpace[i] = new AliL3HistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi); | |
107 | } | |
b2a02bce | 108 | } |
109 | ||
110 | void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t pt_min, | |
111 | Int_t nybin,Float_t phimin,Float_t phimax) | |
e26acabd | 112 | { |
3ef466c5 | 113 | //Create the histograms (parameter space). |
114 | //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis). | |
115 | //The arguments give the range and binning; | |
116 | //nxbin = #bins in kappa | |
117 | //nybin = #bins in phi0 | |
118 | //pt_min = mimium Pt of track (corresponding to maximum kappa) | |
b2a02bce | 119 | //phi_min = mimimum phi0 |
120 | //phi_max = maximum phi0 | |
3ef466c5 | 121 | |
08a66445 | 122 | Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min; |
b2a02bce | 123 | //Double_t torad = AliL3Transform::Pi()/180; |
124 | ||
125 | CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/); | |
e26acabd | 126 | } |
127 | ||
b2a02bce | 128 | void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax, |
129 | Int_t nybin,Float_t ymin,Float_t ymax) | |
4cafa5fc | 130 | { |
4fc9a6a4 | 131 | |
c52cf5d8 | 132 | fParamSpace = new AliL3Histogram*[GetNEtaSegments()]; |
4cafa5fc | 133 | |
99e7186b | 134 | Char_t histname[256]; |
c52cf5d8 | 135 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
4de874d1 | 136 | { |
99e7186b | 137 | sprintf(histname,"paramspace_%d",i); |
b2a02bce | 138 | //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.5,1.5,0.05,nybin,ymin,ymax); |
99e7186b | 139 | fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax); |
4de874d1 | 140 | } |
1eef4efc | 141 | |
142 | #ifdef do_mc | |
c980a70c | 143 | if(fDoMC) |
144 | { | |
88905d4d | 145 | AliL3Histogram *hist = fParamSpace[0]; |
146 | Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2); | |
3e87ef69 | 147 | cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl; |
c980a70c | 148 | fTrackID = new TrackIndex*[GetNEtaSegments()]; |
149 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
150 | fTrackID[i] = new TrackIndex[ncells]; | |
151 | } | |
1eef4efc | 152 | #endif |
4cafa5fc | 153 | } |
154 | ||
e26acabd | 155 | void AliL3HoughTransformer::Reset() |
156 | { | |
48f3c46f | 157 | //Reset all the histograms |
e26acabd | 158 | |
159 | if(!fParamSpace) | |
160 | { | |
161 | LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms") | |
162 | <<"No histograms to reset"<<ENDLOG; | |
163 | return; | |
164 | } | |
165 | ||
c52cf5d8 | 166 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
e26acabd | 167 | fParamSpace[i]->Reset(); |
b2a02bce | 168 | |
1eef4efc | 169 | #ifdef do_mc |
c980a70c | 170 | if(fDoMC) |
171 | { | |
172 | AliL3Histogram *hist = fParamSpace[0]; | |
173 | Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2); | |
174 | for(Int_t i=0; i<GetNEtaSegments(); i++) | |
175 | memset(fTrackID[i],0,ncells*sizeof(TrackIndex)); | |
176 | } | |
1eef4efc | 177 | #endif |
e26acabd | 178 | } |
179 | ||
b5a207b4 | 180 | Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta) |
181 | { | |
3ef466c5 | 182 | //Return the histogram index of the corresponding eta. |
183 | ||
c52cf5d8 | 184 | Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments(); |
185 | Double_t index = (eta-GetEtaMin())/etaslice; | |
b5a207b4 | 186 | return (Int_t)index; |
187 | } | |
188 | ||
3e87ef69 | 189 | void AliL3HoughTransformer::GetEtaIndexes(Double_t eta,Int_t *indexes) |
190 | { | |
191 | //Return histogram indexes in case of overlapping etaslices. | |
192 | ||
193 | Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments(); | |
194 | Int_t index = (Int_t)((eta-GetEtaMin())/etaslice); | |
195 | if(index%2 == 0) | |
196 | { | |
197 | indexes[0] = index; | |
198 | indexes[1] = index - 1; | |
199 | } | |
200 | else | |
201 | { | |
202 | indexes[0] = index - 1; | |
203 | indexes[1] = index; | |
204 | } | |
205 | } | |
206 | ||
7646f3c3 | 207 | inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index) |
208 | { | |
209 | if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0) | |
210 | return 0; | |
211 | if(!fParamSpace[eta_index]) | |
212 | return 0; | |
213 | return fParamSpace[eta_index]; | |
214 | } | |
215 | ||
afd8fed4 | 216 | Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice) |
7646f3c3 | 217 | { |
218 | Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments(); | |
3e87ef69 | 219 | Double_t eta=0; |
220 | if(fEtaOverlap) | |
221 | { | |
222 | Int_t index = eta_index + 1; | |
223 | eta=(Double_t)((index)*eta_slice); | |
224 | } | |
225 | else | |
226 | eta=(Double_t)((eta_index+0.5)*eta_slice); | |
7646f3c3 | 227 | return eta; |
228 | } | |
229 | ||
4fc9a6a4 | 230 | void AliL3HoughTransformer::TransformCircle() |
4de874d1 | 231 | { |
99e7186b | 232 | //Transform the input data with a circle HT. |
3ef466c5 | 233 | //The function loops over all the data, and transforms each pixel with the equations: |
234 | // | |
235 | //kappa = 2/R*sin(phi - phi0) | |
236 | // | |
237 | //where R = sqrt(x*x +y*y), and phi = arctan(y/x) | |
238 | // | |
239 | //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find | |
3e87ef69 | 240 | //which histogram in which the pixel should be transformed, the eta-value is calculated |
3ef466c5 | 241 | //and the proper histogram index is found by GetEtaIndex(eta). |
242 | ||
52a2a604 | 243 | |
c52cf5d8 | 244 | AliL3DigitRowData *tempPt = GetDataPointer(); |
245 | if(!tempPt) | |
4de874d1 | 246 | { |
237d3f5c | 247 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data") |
248 | <<"No input data "<<ENDLOG; | |
249 | return; | |
250 | } | |
99e7186b | 251 | |
3ef466c5 | 252 | //Loop over the padrows: |
26abc209 | 253 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) |
4cafa5fc | 254 | { |
3ef466c5 | 255 | //Get the data on this padrow: |
99e7186b | 256 | AliL3DigitData *digPt = tempPt->fDigitData; |
257 | if(i != (Int_t)tempPt->fRow) | |
4cafa5fc | 258 | { |
d96f6a4a | 259 | cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl; |
99e7186b | 260 | continue; |
261 | } | |
44c7f8de | 262 | |
3ef466c5 | 263 | //Loop over the data on this padrow: |
99e7186b | 264 | for(UInt_t j=0; j<tempPt->fNDigit; j++) |
265 | { | |
266 | UShort_t charge = digPt[j].fCharge; | |
267 | UChar_t pad = digPt[j].fPad; | |
268 | UShort_t time = digPt[j].fTime; | |
b2a02bce | 269 | if((Int_t)charge <= GetLowerThreshold()) |
99e7186b | 270 | continue; |
44c7f8de | 271 | |
b2a02bce | 272 | if((Int_t)charge > GetUpperThreshold()) |
273 | charge = GetUpperThreshold(); | |
274 | ||
4cafa5fc | 275 | Int_t sector,row; |
99e7186b | 276 | Float_t xyz[3]; |
3e87ef69 | 277 | |
3ef466c5 | 278 | //Transform data to local cartesian coordinates: |
4ab9f8f0 | 279 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); |
280 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
3e87ef69 | 281 | |
3ef466c5 | 282 | //Calculate the eta: |
4ab9f8f0 | 283 | Double_t eta = AliL3Transform::GetEta(xyz); |
3ef466c5 | 284 | |
285 | //Get the corresponding index, which determines which histogram to fill: | |
c52cf5d8 | 286 | Int_t eta_index = GetEtaIndex(eta); |
3e87ef69 | 287 | |
c52cf5d8 | 288 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) |
99e7186b | 289 | continue; |
4cafa5fc | 290 | |
e26acabd | 291 | //Get the correct histogrampointer: |
99e7186b | 292 | AliL3Histogram *hist = fParamSpace[eta_index]; |
293 | if(!hist) | |
4cafa5fc | 294 | { |
3e87ef69 | 295 | cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<eta_index<<endl; |
99e7186b | 296 | continue; |
4cafa5fc | 297 | } |
3e87ef69 | 298 | |
3ef466c5 | 299 | //Do the transformation: |
300 | Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
4ab9f8f0 | 301 | Float_t phi = AliL3Transform::GetPhi(xyz); |
4de874d1 | 302 | |
b2a02bce | 303 | |
99e7186b | 304 | //Fill the histogram along the phirange |
305 | for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++) | |
4de874d1 | 306 | { |
99e7186b | 307 | Float_t phi0 = hist->GetBinCenterY(b); |
308 | Float_t kappa = 2*sin(phi - phi0)/R; | |
b2a02bce | 309 | //hist->Fill(kappa,phi0,(int)rint(log((Float_t)charge))); |
3e87ef69 | 310 | hist->Fill(kappa,phi0,charge); |
311 | //hist->Fill(kappa,phi0,1); | |
1eef4efc | 312 | #ifdef do_mc |
c980a70c | 313 | if(fDoMC) |
1eef4efc | 314 | { |
c980a70c | 315 | Int_t bin = hist->FindBin(kappa,phi0); |
316 | for(Int_t t=0; t<3; t++) | |
317 | { | |
318 | Int_t label = digPt[j].fTrackID[t]; | |
319 | if(label < 0) break; | |
320 | UInt_t c; | |
321 | for(c=0; c<MaxTrack; c++) | |
322 | if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0) | |
323 | break; | |
324 | if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl; | |
325 | fTrackID[eta_index][bin].fLabel[c] = label; | |
326 | fTrackID[eta_index][bin].fNHits[c]++; | |
327 | } | |
1eef4efc | 328 | } |
329 | #endif | |
4de874d1 | 330 | } |
4de874d1 | 331 | } |
3ef466c5 | 332 | |
333 | //Move the data pointer to the next padrow: | |
48f3c46f | 334 | AliL3MemHandler::UpdateRowPointer(tempPt); |
4de874d1 | 335 | } |
4de874d1 | 336 | } |
337 | ||
3e87ef69 | 338 | struct Digit { |
339 | Int_t row; | |
340 | Double_t r; | |
341 | Double_t phi; | |
342 | Int_t charge; | |
343 | Digit *next; | |
344 | }; | |
1f1942b8 | 345 | |
3e87ef69 | 346 | struct EtaContainer { |
347 | Digit *first; | |
348 | Digit *last; | |
349 | }; | |
1f1942b8 | 350 | |
3e87ef69 | 351 | void AliL3HoughTransformer::TransformCircleC(Int_t *row_range,Int_t every) |
b1886074 | 352 | { |
353 | //Circle transform, using combinations of every 2 points lying | |
354 | //on different padrows and within the same etaslice. | |
355 | ||
c52cf5d8 | 356 | AliL3DigitRowData *tempPt = GetDataPointer(); |
b1886074 | 357 | if(!tempPt) |
237d3f5c | 358 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data") |
359 | <<"No input data "<<ENDLOG; | |
3e87ef69 | 360 | |
361 | Int_t minrow = AliL3Transform::GetFirstRow(GetPatch()); | |
362 | Int_t maxrow = AliL3Transform::GetLastRow(GetPatch()); | |
363 | if(row_range) | |
364 | { | |
365 | minrow = row_range[0]; | |
366 | maxrow = row_range[1]; | |
367 | if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch())) | |
368 | minrow = AliL3Transform::GetFirstRow(GetPatch()); | |
369 | if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch())) | |
370 | maxrow = AliL3Transform::GetLastRow(GetPatch()); | |
371 | if(minrow > maxrow || maxrow==minrow) | |
372 | { | |
373 | cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl; | |
374 | return; | |
375 | } | |
376 | } | |
377 | else | |
378 | { | |
379 | minrow = AliL3Transform::GetFirstRow(GetPatch()); | |
380 | maxrow = AliL3Transform::GetLastRow(GetPatch()); | |
381 | } | |
382 | ||
3ef466c5 | 383 | Int_t counter=0; |
26abc209 | 384 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) |
b1886074 | 385 | { |
3ef466c5 | 386 | counter += tempPt->fNDigit; |
b1886074 | 387 | AliL3MemHandler::UpdateRowPointer(tempPt); |
388 | } | |
3ef466c5 | 389 | |
3e87ef69 | 390 | Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1); |
391 | EtaContainer *etaPt = new EtaContainer[bound]; | |
392 | memset(etaPt,0,bound*sizeof(EtaContainer)); | |
3ef466c5 | 393 | |
394 | Digit *digits = new Digit[counter]; | |
395 | cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl; | |
3e87ef69 | 396 | memset(digits,0,counter*sizeof(Digit)); |
397 | ||
3ef466c5 | 398 | Int_t sector,row,tot_charge,pad,time,charge; |
b1886074 | 399 | Double_t r1,r2,phi1,phi2,eta,kappa,phi_0; |
b1886074 | 400 | Float_t xyz[3]; |
3ef466c5 | 401 | |
402 | counter=0; | |
c52cf5d8 | 403 | tempPt = GetDataPointer(); |
3ef466c5 | 404 | |
3e87ef69 | 405 | cout<<"Calculating digits in patch "<<GetPatch()<<endl; |
26abc209 | 406 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) |
b1886074 | 407 | { |
3ef466c5 | 408 | AliL3DigitData *digPt = tempPt->fDigitData; |
409 | for(UInt_t di=0; di<tempPt->fNDigit; di++) | |
b1886074 | 410 | { |
3ef466c5 | 411 | charge = digPt[di].fCharge; |
b1886074 | 412 | pad = digPt[di].fPad; |
413 | time = digPt[di].fTime; | |
4ab9f8f0 | 414 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); |
415 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
416 | eta = AliL3Transform::GetEta(xyz); | |
3e87ef69 | 417 | |
3ef466c5 | 418 | digits[counter].row = i; |
419 | digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
420 | digits[counter].phi = atan2(xyz[1],xyz[0]); | |
3ef466c5 | 421 | digits[counter].charge = charge; |
3e87ef69 | 422 | |
423 | if(!fEtaOverlap) | |
c980a70c | 424 | { |
3e87ef69 | 425 | Int_t eta_index = GetEtaIndex(eta); |
426 | ||
427 | Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index; | |
428 | ||
429 | if(index > 0 && index < bound) | |
430 | { | |
431 | if(etaPt[index].first == 0) | |
432 | etaPt[index].first = &digits[counter]; | |
433 | else | |
434 | (etaPt[index].last)->next = &digits[counter]; | |
435 | etaPt[index].last = &digits[counter]; | |
436 | } | |
c980a70c | 437 | } |
3e87ef69 | 438 | else |
439 | { | |
440 | Int_t eta_index[2]; | |
441 | GetEtaIndexes(eta,eta_index); | |
442 | Int_t index[2]; | |
443 | index[0] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index[0]; | |
444 | index[1] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index[1]; | |
445 | if(index[0] == index[1]) | |
446 | { | |
447 | cerr<<"Same etaindexes "<<index[0]<<" "<<index[1]<<endl; | |
448 | exit(5); | |
449 | } | |
450 | ||
451 | Int_t ind = index[0]; | |
452 | if(ind > 0 && ind < bound) | |
453 | { | |
454 | if(etaPt[ind].first == 0) | |
455 | etaPt[ind].first = &digits[counter]; | |
456 | else | |
457 | (etaPt[ind].last)->next = &digits[counter]; | |
458 | etaPt[ind].last = &digits[counter]; | |
459 | } | |
460 | ||
461 | ind = index[1]; | |
462 | if(ind > 0 && ind < bound) | |
463 | { | |
464 | if(etaPt[ind].first == 0) | |
465 | etaPt[ind].first = &digits[counter]; | |
466 | else | |
467 | (etaPt[ind].last)->next = &digits[counter]; | |
468 | etaPt[ind].last = &digits[counter]; | |
469 | } | |
470 | } | |
471 | ||
3ef466c5 | 472 | counter++; |
473 | } | |
474 | AliL3MemHandler::UpdateRowPointer(tempPt); | |
475 | } | |
476 | ||
3e87ef69 | 477 | cout<<"Doing the combinatorics"<<endl; |
478 | ||
479 | Digit *dPt1,*dPt2; | |
480 | ||
481 | for(Int_t e=0; e<GetNEtaSegments(); e++) | |
3ef466c5 | 482 | { |
3e87ef69 | 483 | for(Int_t i=minrow; i<=maxrow; i+=every) |
3ef466c5 | 484 | { |
3e87ef69 | 485 | Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e; |
3ef466c5 | 486 | |
3e87ef69 | 487 | for(dPt1 = (Digit*)etaPt[index1].first; dPt1 != 0; dPt1 = (Digit*)dPt1->next) |
1eef4efc | 488 | { |
3e87ef69 | 489 | for(Int_t j=i+every; j<=maxrow; j+=every) |
1eef4efc | 490 | { |
3e87ef69 | 491 | Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e; |
492 | ||
493 | for(dPt2 = (Digit*)etaPt[index2].first; dPt2 != 0; dPt2 = (Digit*)dPt2->next) | |
1eef4efc | 494 | { |
3e87ef69 | 495 | if(dPt1->row == dPt2->row) |
c980a70c | 496 | { |
3e87ef69 | 497 | cerr<<"same row; indexes "<<index1<<" "<<index2<<endl; |
498 | exit(5); | |
c980a70c | 499 | } |
3e87ef69 | 500 | |
501 | //Get the correct histogrampointer: | |
502 | AliL3Histogram *hist = fParamSpace[e]; | |
503 | if(!hist) | |
504 | { | |
505 | printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i); | |
506 | continue; | |
507 | } | |
508 | ||
509 | //Do the transform: | |
510 | r1 = dPt1->r; | |
511 | phi1 = dPt1->phi; | |
512 | r2 = dPt2->r; | |
513 | phi2 = dPt2->phi; | |
514 | phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) ); | |
515 | kappa = 2*sin(phi2-phi_0)/r2; | |
516 | tot_charge = dPt1->charge + dPt2->charge; | |
517 | hist->Fill(kappa,phi_0,tot_charge); | |
518 | ||
1eef4efc | 519 | } |
520 | } | |
521 | } | |
b1886074 | 522 | } |
523 | } | |
3e87ef69 | 524 | |
525 | cout<<"done"<<endl; | |
526 | delete [] etaPt; | |
3ef466c5 | 527 | delete [] digits; |
3e87ef69 | 528 | |
b1886074 | 529 | } |
530 | ||
3e87ef69 | 531 | void AliL3HoughTransformer::TransformLine(Int_t *row_range,Float_t *phirange) |
4fc9a6a4 | 532 | { |
533 | //Do a line transform on the data. | |
534 | ||
3e87ef69 | 535 | |
c52cf5d8 | 536 | AliL3DigitRowData *tempPt = GetDataPointer(); |
537 | if(!tempPt) | |
4fc9a6a4 | 538 | { |
237d3f5c | 539 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data") |
540 | <<"No input data "<<ENDLOG; | |
541 | return; | |
542 | } | |
3e87ef69 | 543 | |
544 | Int_t minrow = AliL3Transform::GetFirstRow(GetPatch()); | |
545 | Int_t maxrow = AliL3Transform::GetLastRow(GetPatch()); | |
546 | if(row_range) | |
547 | { | |
548 | minrow = row_range[0]; | |
549 | maxrow = row_range[1]; | |
550 | if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch())) | |
551 | minrow = AliL3Transform::GetFirstRow(GetPatch()); | |
552 | if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch())) | |
553 | maxrow = AliL3Transform::GetLastRow(GetPatch()); | |
554 | if(minrow > maxrow || maxrow==minrow) | |
555 | { | |
556 | cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl; | |
557 | return; | |
558 | } | |
559 | } | |
560 | ||
561 | for(Int_t i=minrow; i<=maxrow; i++) | |
4fc9a6a4 | 562 | { |
563 | AliL3DigitData *digPt = tempPt->fDigitData; | |
564 | if(i != (Int_t)tempPt->fRow) | |
565 | { | |
566 | printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n"); | |
567 | continue; | |
568 | } | |
569 | for(UInt_t j=0; j<tempPt->fNDigit; j++) | |
570 | { | |
571 | UShort_t charge = digPt[j].fCharge; | |
572 | UChar_t pad = digPt[j].fPad; | |
573 | UShort_t time = digPt[j].fTime; | |
3bb06991 | 574 | if(charge < GetLowerThreshold()) |
4fc9a6a4 | 575 | continue; |
576 | Int_t sector,row; | |
577 | Float_t xyz[3]; | |
4ab9f8f0 | 578 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); |
579 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
3e87ef69 | 580 | |
581 | if(phirange) | |
582 | { | |
583 | Float_t phi = AliL3Transform::GetPhi(xyz); | |
584 | if(phi < phirange[0] || phi > phirange[1]) | |
585 | continue; | |
586 | } | |
4ab9f8f0 | 587 | Float_t eta = AliL3Transform::GetEta(xyz); |
b5a207b4 | 588 | Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice); |
c52cf5d8 | 589 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) |
4fc9a6a4 | 590 | continue; |
591 | ||
3e87ef69 | 592 | xyz[0] = xyz[0] - AliL3Transform::Row2X(minrow); |
593 | ||
4fc9a6a4 | 594 | //Get the correct histogram: |
595 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
596 | if(!hist) | |
597 | { | |
598 | printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index); | |
599 | continue; | |
600 | } | |
601 | for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++) | |
602 | { | |
603 | Double_t theta = hist->GetBinCenterX(xbin); | |
604 | Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta); | |
605 | hist->Fill(theta,rho,charge); | |
606 | } | |
607 | } | |
48f3c46f | 608 | AliL3MemHandler::UpdateRowPointer(tempPt); |
4fc9a6a4 | 609 | } |
610 | ||
611 | } | |
1eef4efc | 612 | |
3e87ef69 | 613 | struct LDigit { |
614 | Int_t row; | |
615 | Int_t charge; | |
616 | Float_t y; | |
617 | LDigit *next; | |
618 | }; | |
619 | struct LEtaContainer { | |
620 | LDigit *first; | |
621 | LDigit *last; | |
622 | }; | |
623 | void AliL3HoughTransformer::TransformLineC(Int_t *rowrange,Float_t *phirange) | |
624 | { | |
625 | AliL3DigitRowData *tempPt = GetDataPointer(); | |
626 | if(!tempPt) | |
627 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data") | |
628 | <<"No input data "<<ENDLOG; | |
629 | ||
630 | ||
631 | Int_t counter=0; | |
632 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) | |
633 | { | |
634 | counter += tempPt->fNDigit; | |
635 | AliL3MemHandler::UpdateRowPointer(tempPt); | |
636 | } | |
637 | ||
638 | Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1); | |
639 | LEtaContainer *etaPt = new LEtaContainer[bound]; | |
640 | memset(etaPt,0,bound*sizeof(LEtaContainer)); | |
641 | ||
642 | LDigit *digits = new LDigit[counter]; | |
643 | cout<<"Allocating "<<counter*sizeof(LDigit)<<" bytes to digitsarray"<<endl; | |
644 | memset(digits,0,counter*sizeof(LDigit)); | |
645 | ||
646 | Int_t sector,row; | |
647 | Float_t xyz[3]; | |
648 | ||
649 | counter=0; | |
650 | tempPt = GetDataPointer(); | |
651 | ||
652 | cout<<"Calculating digits in patch "<<GetPatch()<<endl; | |
653 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) | |
654 | { | |
655 | AliL3DigitData *digPt = tempPt->fDigitData; | |
656 | for(UInt_t di=0; di<tempPt->fNDigit; di++) | |
657 | { | |
658 | Int_t charge = digPt[di].fCharge; | |
659 | Int_t pad = digPt[di].fPad; | |
660 | Int_t time = digPt[di].fTime; | |
661 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); | |
662 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
663 | Double_t eta = AliL3Transform::GetEta(xyz); | |
664 | ||
665 | Float_t phi = atan2(xyz[1],xyz[0]); | |
666 | if(phi < phirange[0] || phi > phirange[1]) continue; | |
667 | ||
668 | digits[counter].row = i; | |
669 | digits[counter].y = xyz[1]; | |
670 | digits[counter].charge = charge; | |
671 | ||
672 | Int_t eta_index = GetEtaIndex(eta); | |
673 | Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index; | |
674 | ||
675 | if(index > 0 && index < bound) | |
676 | { | |
677 | if(etaPt[index].first == 0) | |
678 | etaPt[index].first = &digits[counter]; | |
679 | else | |
680 | (etaPt[index].last)->next = &digits[counter]; | |
681 | etaPt[index].last = &digits[counter]; | |
682 | } | |
683 | counter++; | |
684 | } | |
685 | AliL3MemHandler::UpdateRowPointer(tempPt); | |
686 | } | |
687 | ||
688 | cout<<"Doing the combinatorics"<<endl; | |
689 | ||
690 | LDigit *dPt1,*dPt2; | |
691 | ||
692 | for(Int_t e=0; e<GetNEtaSegments(); e++) | |
693 | { | |
694 | for(Int_t i=rowrange[0]; i<=rowrange[1]; i++) | |
695 | { | |
696 | Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e; | |
697 | ||
698 | for(dPt1 = (LDigit*)etaPt[index1].first; dPt1 != 0; dPt1 = (LDigit*)dPt1->next) | |
699 | { | |
700 | for(Int_t j=i+1; j<=rowrange[1]; j++) | |
701 | { | |
702 | Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e; | |
703 | ||
704 | for(dPt2 = (LDigit*)etaPt[index2].first; dPt2 != 0; dPt2 = (LDigit*)dPt2->next) | |
705 | { | |
706 | if(dPt1->row == dPt2->row) | |
707 | { | |
708 | cerr<<"same row; indexes "<<index1<<" "<<index2<<endl; | |
709 | exit(5); | |
710 | } | |
711 | ||
712 | //Get the correct histogrampointer: | |
713 | AliL3Histogram *hist = fParamSpace[e]; | |
714 | if(!hist) | |
715 | { | |
716 | printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i); | |
717 | continue; | |
718 | } | |
719 | ||
720 | //Do the transform: | |
721 | float x1 = AliL3Transform::Row2X(dPt1->row) - AliL3Transform::Row2X(rowrange[0]); | |
722 | float x2 = AliL3Transform::Row2X(dPt2->row) - AliL3Transform::Row2X(rowrange[0]); | |
723 | float y1 = dPt1->y; | |
724 | float y2 = dPt2->y; | |
725 | float theta = atan2(x2-x1,y1-y2); | |
726 | float rho = x1*cos(theta)+y1*sin(theta); | |
727 | hist->Fill(theta,rho,1);//dPt1->charge+dPt2->charge); | |
728 | } | |
729 | } | |
730 | } | |
731 | } | |
732 | } | |
733 | ||
734 | cout<<"done"<<endl; | |
735 | delete [] etaPt; | |
736 | delete [] digits; | |
737 | } | |
738 | ||
1eef4efc | 739 | Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi) |
740 | { | |
c980a70c | 741 | if(!fDoMC) |
742 | { | |
743 | cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl; | |
744 | return -1; | |
745 | } | |
746 | ||
1eef4efc | 747 | #ifdef do_mc |
748 | if(eta_index < 0 || eta_index > GetNEtaSegments()) | |
749 | { | |
750 | cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl; | |
751 | return -1; | |
752 | } | |
753 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
754 | Int_t bin = hist->FindBin(kappa,psi); | |
755 | Int_t label=-1; | |
756 | Int_t max=0; | |
757 | for(UInt_t i=0; i<MaxTrack; i++) | |
758 | { | |
759 | Int_t nhits=fTrackID[eta_index][bin].fNHits[i]; | |
760 | if(nhits == 0) break; | |
761 | if(nhits > max) | |
762 | { | |
763 | max = nhits; | |
764 | label = fTrackID[eta_index][bin].fLabel[i]; | |
765 | } | |
766 | } | |
3e87ef69 | 767 | //nhits = max; |
1eef4efc | 768 | return label; |
769 | #endif | |
770 | cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl; | |
771 | return -1; | |
772 | } | |
773 |