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