]>
Commit | Line | Data |
---|---|---|
b5a207b4 | 1 | //$Id$ |
2 | ||
b1886074 | 3 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> |
4 | //*-- Copyright © ASV | |
4cafa5fc | 5 | |
48f3c46f | 6 | #include "AliL3MemHandler.h" |
e26acabd | 7 | #include "AliL3Logging.h" |
99e7186b | 8 | #include "AliL3HoughTransformer.h" |
4de874d1 | 9 | #include "AliL3Defs.h" |
10 | #include "AliL3Transform.h" | |
4cafa5fc | 11 | #include "AliL3DigitData.h" |
99e7186b | 12 | #include "AliL3Histogram.h" |
4de874d1 | 13 | |
b1886074 | 14 | //_____________________________________________________________ |
15 | // AliL3HoughTransformer | |
16 | // | |
17 | // Hough transformation class | |
3ef466c5 | 18 | // |
b1886074 | 19 | |
4de874d1 | 20 | ClassImp(AliL3HoughTransformer) |
21 | ||
22 | AliL3HoughTransformer::AliL3HoughTransformer() | |
23 | { | |
24 | //Default constructor | |
c52cf5d8 | 25 | fParamSpace = 0; |
4de874d1 | 26 | } |
27 | ||
c52cf5d8 | 28 | AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments) |
52a2a604 | 29 | { |
3ef466c5 | 30 | //Normal constructor |
237d3f5c | 31 | fParamSpace = 0; |
52a2a604 | 32 | } |
33 | ||
4de874d1 | 34 | AliL3HoughTransformer::~AliL3HoughTransformer() |
35 | { | |
99e7186b | 36 | DeleteHistograms(); |
4cafa5fc | 37 | } |
38 | ||
99e7186b | 39 | void AliL3HoughTransformer::DeleteHistograms() |
4cafa5fc | 40 | { |
99e7186b | 41 | if(!fParamSpace) |
42 | return; | |
c52cf5d8 | 43 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
4fc9a6a4 | 44 | { |
45 | if(!fParamSpace[i]) continue; | |
46 | delete fParamSpace[i]; | |
47 | } | |
99e7186b | 48 | delete [] fParamSpace; |
4cafa5fc | 49 | } |
50 | ||
e26acabd | 51 | void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min, |
52 | Int_t nybin,Double_t phimin,Double_t phimax) | |
53 | { | |
3ef466c5 | 54 | //Create the histograms (parameter space). |
55 | //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis). | |
56 | //The arguments give the range and binning; | |
57 | //nxbin = #bins in kappa | |
58 | //nybin = #bins in phi0 | |
59 | //pt_min = mimium Pt of track (corresponding to maximum kappa) | |
60 | //phi_min = mimimum phi0 (degrees) | |
61 | //phi_max = maximum phi0 (degrees) | |
62 | ||
e26acabd | 63 | Double_t bfact = 0.0029980; |
64 | Double_t bfield = 0.2; | |
65 | Double_t x = bfact*bfield/pt_min; | |
b5a207b4 | 66 | CreateHistograms(nxbin,-1.*x,x,nybin,phimin*ToRad,phimax*ToRad); |
e26acabd | 67 | } |
68 | ||
99e7186b | 69 | void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax, |
70 | Int_t nybin,Double_t ymin,Double_t ymax) | |
4cafa5fc | 71 | { |
4fc9a6a4 | 72 | |
c52cf5d8 | 73 | fParamSpace = new AliL3Histogram*[GetNEtaSegments()]; |
4cafa5fc | 74 | |
99e7186b | 75 | Char_t histname[256]; |
c52cf5d8 | 76 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
4de874d1 | 77 | { |
99e7186b | 78 | sprintf(histname,"paramspace_%d",i); |
79 | fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax); | |
4de874d1 | 80 | } |
4cafa5fc | 81 | } |
82 | ||
e26acabd | 83 | void AliL3HoughTransformer::Reset() |
84 | { | |
48f3c46f | 85 | //Reset all the histograms |
e26acabd | 86 | |
87 | if(!fParamSpace) | |
88 | { | |
89 | LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms") | |
90 | <<"No histograms to reset"<<ENDLOG; | |
91 | return; | |
92 | } | |
93 | ||
c52cf5d8 | 94 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
e26acabd | 95 | fParamSpace[i]->Reset(); |
e26acabd | 96 | } |
97 | ||
4de874d1 | 98 | |
b5a207b4 | 99 | Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta) |
100 | { | |
3ef466c5 | 101 | //Return the histogram index of the corresponding eta. |
102 | ||
c52cf5d8 | 103 | Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments(); |
104 | Double_t index = (eta-GetEtaMin())/etaslice; | |
b5a207b4 | 105 | return (Int_t)index; |
106 | } | |
107 | ||
4fc9a6a4 | 108 | void AliL3HoughTransformer::TransformCircle() |
4de874d1 | 109 | { |
99e7186b | 110 | //Transform the input data with a circle HT. |
3ef466c5 | 111 | //The function loops over all the data, and transforms each pixel with the equations: |
112 | // | |
113 | //kappa = 2/R*sin(phi - phi0) | |
114 | // | |
115 | //where R = sqrt(x*x +y*y), and phi = arctan(y/x) | |
116 | // | |
117 | //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find | |
118 | //which histogram in which the pixel should be transformed, the eta-value is calcluated | |
119 | //and the proper histogram index is found by GetEtaIndex(eta). | |
120 | ||
52a2a604 | 121 | |
c52cf5d8 | 122 | AliL3DigitRowData *tempPt = GetDataPointer(); |
123 | if(!tempPt) | |
4de874d1 | 124 | { |
237d3f5c | 125 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data") |
126 | <<"No input data "<<ENDLOG; | |
127 | return; | |
128 | } | |
129 | if(!fTransform) | |
130 | { | |
131 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Transformer") | |
132 | <<"No AliL3Transform object"<<ENDLOG; | |
99e7186b | 133 | return; |
4de874d1 | 134 | } |
99e7186b | 135 | |
3ef466c5 | 136 | //Loop over the padrows: |
c52cf5d8 | 137 | for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++) |
4cafa5fc | 138 | { |
3ef466c5 | 139 | //Get the data on this padrow: |
99e7186b | 140 | AliL3DigitData *digPt = tempPt->fDigitData; |
141 | if(i != (Int_t)tempPt->fRow) | |
4cafa5fc | 142 | { |
4fc9a6a4 | 143 | printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n"); |
99e7186b | 144 | continue; |
145 | } | |
3ef466c5 | 146 | |
147 | //Loop over the data on this padrow: | |
99e7186b | 148 | for(UInt_t j=0; j<tempPt->fNDigit; j++) |
149 | { | |
150 | UShort_t charge = digPt[j].fCharge; | |
151 | UChar_t pad = digPt[j].fPad; | |
152 | UShort_t time = digPt[j].fTime; | |
c52cf5d8 | 153 | if(charge <= GetThreshold()) |
99e7186b | 154 | continue; |
4cafa5fc | 155 | Int_t sector,row; |
99e7186b | 156 | Float_t xyz[3]; |
3ef466c5 | 157 | |
158 | //Transform data to local cartesian coordinates: | |
c52cf5d8 | 159 | fTransform->Slice2Sector(GetSlice(),i,sector,row); |
99e7186b | 160 | fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); |
3ef466c5 | 161 | |
162 | //Calculate the eta: | |
4fc9a6a4 | 163 | Double_t eta = fTransform->GetEta(xyz); |
3ef466c5 | 164 | |
165 | //Get the corresponding index, which determines which histogram to fill: | |
c52cf5d8 | 166 | Int_t eta_index = GetEtaIndex(eta); |
167 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) | |
99e7186b | 168 | continue; |
4cafa5fc | 169 | |
e26acabd | 170 | //Get the correct histogrampointer: |
99e7186b | 171 | AliL3Histogram *hist = fParamSpace[eta_index]; |
172 | if(!hist) | |
4cafa5fc | 173 | { |
4fc9a6a4 | 174 | printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index); |
99e7186b | 175 | continue; |
4cafa5fc | 176 | } |
4cafa5fc | 177 | |
3ef466c5 | 178 | //Do the transformation: |
179 | Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
99e7186b | 180 | Float_t phi = fTransform->GetPhi(xyz); |
4de874d1 | 181 | |
99e7186b | 182 | //Fill the histogram along the phirange |
183 | for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++) | |
4de874d1 | 184 | { |
99e7186b | 185 | Float_t phi0 = hist->GetBinCenterY(b); |
186 | Float_t kappa = 2*sin(phi - phi0)/R; | |
187 | hist->Fill(kappa,phi0,charge); | |
4de874d1 | 188 | } |
4de874d1 | 189 | } |
3ef466c5 | 190 | |
191 | //Move the data pointer to the next padrow: | |
48f3c46f | 192 | AliL3MemHandler::UpdateRowPointer(tempPt); |
4de874d1 | 193 | } |
4de874d1 | 194 | } |
195 | ||
3ef466c5 | 196 | void AliL3HoughTransformer::TransformCircleC(Int_t row_range) |
b1886074 | 197 | { |
198 | //Circle transform, using combinations of every 2 points lying | |
199 | //on different padrows and within the same etaslice. | |
200 | ||
c52cf5d8 | 201 | AliL3DigitRowData *tempPt = GetDataPointer(); |
b1886074 | 202 | if(!tempPt) |
237d3f5c | 203 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data") |
204 | <<"No input data "<<ENDLOG; | |
205 | ||
206 | if(!fTransform) | |
207 | { | |
208 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Transformer") | |
209 | <<"No AliL3Transform object"<<ENDLOG; | |
210 | return; | |
211 | } | |
b1886074 | 212 | |
3ef466c5 | 213 | Int_t counter=0; |
c52cf5d8 | 214 | for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++) |
b1886074 | 215 | { |
3ef466c5 | 216 | counter += tempPt->fNDigit; |
b1886074 | 217 | AliL3MemHandler::UpdateRowPointer(tempPt); |
218 | } | |
3ef466c5 | 219 | |
220 | struct Digit { | |
221 | Int_t row; | |
222 | Double_t r; | |
223 | Double_t phi; | |
224 | Int_t eta_index; | |
225 | Int_t charge; | |
226 | }; | |
227 | ||
228 | Digit *digits = new Digit[counter]; | |
229 | cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl; | |
230 | ||
231 | Int_t total_digits=counter; | |
232 | Int_t sector,row,tot_charge,pad,time,charge; | |
b1886074 | 233 | Double_t r1,r2,phi1,phi2,eta,kappa,phi_0; |
b1886074 | 234 | Float_t xyz[3]; |
3ef466c5 | 235 | |
236 | counter=0; | |
c52cf5d8 | 237 | tempPt = GetDataPointer(); |
3ef466c5 | 238 | |
c52cf5d8 | 239 | for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++) |
b1886074 | 240 | { |
3ef466c5 | 241 | AliL3DigitData *digPt = tempPt->fDigitData; |
242 | for(UInt_t di=0; di<tempPt->fNDigit; di++) | |
b1886074 | 243 | { |
3ef466c5 | 244 | charge = digPt[di].fCharge; |
b1886074 | 245 | pad = digPt[di].fPad; |
246 | time = digPt[di].fTime; | |
c52cf5d8 | 247 | fTransform->Slice2Sector(GetSlice(),i,sector,row); |
b1886074 | 248 | fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); |
249 | eta = fTransform->GetEta(xyz); | |
3ef466c5 | 250 | digits[counter].row = i; |
251 | digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
252 | digits[counter].phi = atan2(xyz[1],xyz[0]); | |
253 | digits[counter].eta_index = GetEtaIndex(eta); | |
254 | digits[counter].charge = charge; | |
255 | counter++; | |
256 | } | |
257 | AliL3MemHandler::UpdateRowPointer(tempPt); | |
258 | } | |
259 | ||
260 | for(Int_t i=0; i<total_digits; i++) | |
261 | { | |
c52cf5d8 | 262 | if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue; |
3ef466c5 | 263 | Int_t ind = digits[i].eta_index; |
264 | ||
265 | for(Int_t j=i+1; j<total_digits; j++) | |
266 | { | |
267 | if(digits[i].row == digits[j].row) continue; | |
268 | if(digits[i].eta_index != digits[j].eta_index) continue; | |
269 | if(digits[i].row + row_range < digits[j].row) break; | |
b1886074 | 270 | |
271 | //Get the correct histogrampointer: | |
3ef466c5 | 272 | AliL3Histogram *hist = fParamSpace[ind]; |
b1886074 | 273 | if(!hist) |
274 | { | |
3ef466c5 | 275 | printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind); |
b1886074 | 276 | continue; |
277 | } | |
3ef466c5 | 278 | |
279 | r1 = digits[i].r; | |
280 | phi1 = digits[i].phi; | |
281 | r2 = digits[j].r; | |
282 | phi2 = digits[j].phi; | |
283 | phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) ); | |
284 | kappa = 2*sin(phi2-phi_0)/r2; | |
285 | tot_charge = digits[i].charge + digits[j].charge; | |
286 | hist->Fill(kappa,phi_0,tot_charge); | |
b1886074 | 287 | } |
288 | } | |
3ef466c5 | 289 | delete [] digits; |
b1886074 | 290 | } |
291 | ||
4fc9a6a4 | 292 | void AliL3HoughTransformer::TransformLine() |
293 | { | |
294 | //Do a line transform on the data. | |
295 | ||
296 | ||
c52cf5d8 | 297 | AliL3DigitRowData *tempPt = GetDataPointer(); |
298 | if(!tempPt) | |
4fc9a6a4 | 299 | { |
237d3f5c | 300 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data") |
301 | <<"No input data "<<ENDLOG; | |
302 | return; | |
303 | } | |
304 | if(!fTransform) | |
305 | { | |
306 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Transformer") | |
307 | <<"No AliL3Transform object"<<ENDLOG; | |
4fc9a6a4 | 308 | return; |
309 | } | |
310 | ||
c52cf5d8 | 311 | for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++) |
4fc9a6a4 | 312 | { |
313 | AliL3DigitData *digPt = tempPt->fDigitData; | |
314 | if(i != (Int_t)tempPt->fRow) | |
315 | { | |
316 | printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n"); | |
317 | continue; | |
318 | } | |
319 | for(UInt_t j=0; j<tempPt->fNDigit; j++) | |
320 | { | |
321 | UShort_t charge = digPt[j].fCharge; | |
322 | UChar_t pad = digPt[j].fPad; | |
323 | UShort_t time = digPt[j].fTime; | |
c52cf5d8 | 324 | if(charge < GetThreshold()) |
4fc9a6a4 | 325 | continue; |
326 | Int_t sector,row; | |
327 | Float_t xyz[3]; | |
c52cf5d8 | 328 | fTransform->Slice2Sector(GetSlice(),i,sector,row); |
4fc9a6a4 | 329 | fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); |
330 | Float_t eta = fTransform->GetEta(xyz); | |
b5a207b4 | 331 | Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice); |
c52cf5d8 | 332 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) |
4fc9a6a4 | 333 | continue; |
334 | ||
335 | //Get the correct histogram: | |
336 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
337 | if(!hist) | |
338 | { | |
339 | printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index); | |
340 | continue; | |
341 | } | |
342 | for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++) | |
343 | { | |
344 | Double_t theta = hist->GetBinCenterX(xbin); | |
345 | Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta); | |
346 | hist->Fill(theta,rho,charge); | |
347 | } | |
348 | } | |
48f3c46f | 349 | AliL3MemHandler::UpdateRowPointer(tempPt); |
4fc9a6a4 | 350 | } |
351 | ||
352 | } |