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