]>
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 "AliL3Transform.h" |
4cafa5fc | 10 | #include "AliL3DigitData.h" |
99e7186b | 11 | #include "AliL3Histogram.h" |
4de874d1 | 12 | |
b1886074 | 13 | //_____________________________________________________________ |
14 | // AliL3HoughTransformer | |
15 | // | |
16 | // Hough transformation class | |
3ef466c5 | 17 | // |
b1886074 | 18 | |
4de874d1 | 19 | ClassImp(AliL3HoughTransformer) |
20 | ||
21 | AliL3HoughTransformer::AliL3HoughTransformer() | |
22 | { | |
23 | //Default constructor | |
c52cf5d8 | 24 | fParamSpace = 0; |
4de874d1 | 25 | } |
26 | ||
c52cf5d8 | 27 | AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments) |
52a2a604 | 28 | { |
3ef466c5 | 29 | //Normal constructor |
237d3f5c | 30 | fParamSpace = 0; |
52a2a604 | 31 | } |
32 | ||
4de874d1 | 33 | AliL3HoughTransformer::~AliL3HoughTransformer() |
34 | { | |
99e7186b | 35 | DeleteHistograms(); |
4cafa5fc | 36 | } |
37 | ||
99e7186b | 38 | void AliL3HoughTransformer::DeleteHistograms() |
4cafa5fc | 39 | { |
99e7186b | 40 | if(!fParamSpace) |
41 | return; | |
c52cf5d8 | 42 | for(Int_t i=0; i<GetNEtaSegments(); i++) |
4fc9a6a4 | 43 | { |
44 | if(!fParamSpace[i]) continue; | |
45 | delete fParamSpace[i]; | |
46 | } | |
99e7186b | 47 | delete [] fParamSpace; |
4cafa5fc | 48 | } |
49 | ||
e26acabd | 50 | void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min, |
51 | Int_t nybin,Double_t phimin,Double_t phimax) | |
52 | { | |
3ef466c5 | 53 | //Create the histograms (parameter space). |
54 | //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis). | |
55 | //The arguments give the range and binning; | |
56 | //nxbin = #bins in kappa | |
57 | //nybin = #bins in phi0 | |
58 | //pt_min = mimium Pt of track (corresponding to maximum kappa) | |
59 | //phi_min = mimimum phi0 (degrees) | |
60 | //phi_max = maximum phi0 (degrees) | |
61 | ||
e26acabd | 62 | Double_t bfact = 0.0029980; |
63 | Double_t bfield = 0.2; | |
64 | Double_t x = bfact*bfield/pt_min; | |
26abc209 | 65 | Double_t torad = AliL3Transform::Pi()/180; |
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 | } | |
99e7186b | 129 | |
3ef466c5 | 130 | //Loop over the padrows: |
26abc209 | 131 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) |
4cafa5fc | 132 | { |
3ef466c5 | 133 | //Get the data on this padrow: |
99e7186b | 134 | AliL3DigitData *digPt = tempPt->fDigitData; |
135 | if(i != (Int_t)tempPt->fRow) | |
4cafa5fc | 136 | { |
4fc9a6a4 | 137 | printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n"); |
99e7186b | 138 | continue; |
139 | } | |
3ef466c5 | 140 | |
141 | //Loop over the data on this padrow: | |
99e7186b | 142 | for(UInt_t j=0; j<tempPt->fNDigit; j++) |
143 | { | |
144 | UShort_t charge = digPt[j].fCharge; | |
145 | UChar_t pad = digPt[j].fPad; | |
146 | UShort_t time = digPt[j].fTime; | |
c52cf5d8 | 147 | if(charge <= GetThreshold()) |
99e7186b | 148 | continue; |
4cafa5fc | 149 | Int_t sector,row; |
99e7186b | 150 | Float_t xyz[3]; |
3ef466c5 | 151 | |
152 | //Transform data to local cartesian coordinates: | |
4ab9f8f0 | 153 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); |
154 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
3ef466c5 | 155 | |
156 | //Calculate the eta: | |
4ab9f8f0 | 157 | Double_t eta = AliL3Transform::GetEta(xyz); |
3ef466c5 | 158 | |
159 | //Get the corresponding index, which determines which histogram to fill: | |
c52cf5d8 | 160 | Int_t eta_index = GetEtaIndex(eta); |
161 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) | |
99e7186b | 162 | continue; |
4cafa5fc | 163 | |
e26acabd | 164 | //Get the correct histogrampointer: |
99e7186b | 165 | AliL3Histogram *hist = fParamSpace[eta_index]; |
166 | if(!hist) | |
4cafa5fc | 167 | { |
4fc9a6a4 | 168 | printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index); |
99e7186b | 169 | continue; |
4cafa5fc | 170 | } |
4cafa5fc | 171 | |
3ef466c5 | 172 | //Do the transformation: |
173 | Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
4ab9f8f0 | 174 | Float_t phi = AliL3Transform::GetPhi(xyz); |
4de874d1 | 175 | |
99e7186b | 176 | //Fill the histogram along the phirange |
177 | for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++) | |
4de874d1 | 178 | { |
99e7186b | 179 | Float_t phi0 = hist->GetBinCenterY(b); |
180 | Float_t kappa = 2*sin(phi - phi0)/R; | |
181 | hist->Fill(kappa,phi0,charge); | |
4de874d1 | 182 | } |
4de874d1 | 183 | } |
3ef466c5 | 184 | |
185 | //Move the data pointer to the next padrow: | |
48f3c46f | 186 | AliL3MemHandler::UpdateRowPointer(tempPt); |
4de874d1 | 187 | } |
4de874d1 | 188 | } |
189 | ||
3ef466c5 | 190 | void AliL3HoughTransformer::TransformCircleC(Int_t row_range) |
b1886074 | 191 | { |
192 | //Circle transform, using combinations of every 2 points lying | |
193 | //on different padrows and within the same etaslice. | |
194 | ||
c52cf5d8 | 195 | AliL3DigitRowData *tempPt = GetDataPointer(); |
b1886074 | 196 | if(!tempPt) |
237d3f5c | 197 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data") |
198 | <<"No input data "<<ENDLOG; | |
199 | ||
3ef466c5 | 200 | Int_t counter=0; |
26abc209 | 201 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) |
b1886074 | 202 | { |
3ef466c5 | 203 | counter += tempPt->fNDigit; |
b1886074 | 204 | AliL3MemHandler::UpdateRowPointer(tempPt); |
205 | } | |
3ef466c5 | 206 | |
207 | struct Digit { | |
208 | Int_t row; | |
209 | Double_t r; | |
210 | Double_t phi; | |
211 | Int_t eta_index; | |
212 | Int_t charge; | |
213 | }; | |
214 | ||
215 | Digit *digits = new Digit[counter]; | |
216 | cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl; | |
217 | ||
218 | Int_t total_digits=counter; | |
219 | Int_t sector,row,tot_charge,pad,time,charge; | |
b1886074 | 220 | Double_t r1,r2,phi1,phi2,eta,kappa,phi_0; |
b1886074 | 221 | Float_t xyz[3]; |
3ef466c5 | 222 | |
223 | counter=0; | |
c52cf5d8 | 224 | tempPt = GetDataPointer(); |
3ef466c5 | 225 | |
26abc209 | 226 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) |
b1886074 | 227 | { |
3ef466c5 | 228 | AliL3DigitData *digPt = tempPt->fDigitData; |
229 | for(UInt_t di=0; di<tempPt->fNDigit; di++) | |
b1886074 | 230 | { |
3ef466c5 | 231 | charge = digPt[di].fCharge; |
b1886074 | 232 | pad = digPt[di].fPad; |
233 | time = digPt[di].fTime; | |
4ab9f8f0 | 234 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); |
235 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
236 | eta = AliL3Transform::GetEta(xyz); | |
3ef466c5 | 237 | digits[counter].row = i; |
238 | digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); | |
239 | digits[counter].phi = atan2(xyz[1],xyz[0]); | |
240 | digits[counter].eta_index = GetEtaIndex(eta); | |
241 | digits[counter].charge = charge; | |
242 | counter++; | |
243 | } | |
244 | AliL3MemHandler::UpdateRowPointer(tempPt); | |
245 | } | |
246 | ||
247 | for(Int_t i=0; i<total_digits; i++) | |
248 | { | |
c52cf5d8 | 249 | if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue; |
3ef466c5 | 250 | Int_t ind = digits[i].eta_index; |
251 | ||
252 | for(Int_t j=i+1; j<total_digits; j++) | |
253 | { | |
254 | if(digits[i].row == digits[j].row) continue; | |
255 | if(digits[i].eta_index != digits[j].eta_index) continue; | |
256 | if(digits[i].row + row_range < digits[j].row) break; | |
b1886074 | 257 | |
258 | //Get the correct histogrampointer: | |
3ef466c5 | 259 | AliL3Histogram *hist = fParamSpace[ind]; |
b1886074 | 260 | if(!hist) |
261 | { | |
3ef466c5 | 262 | printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind); |
b1886074 | 263 | continue; |
264 | } | |
3ef466c5 | 265 | |
266 | r1 = digits[i].r; | |
267 | phi1 = digits[i].phi; | |
268 | r2 = digits[j].r; | |
269 | phi2 = digits[j].phi; | |
270 | phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) ); | |
271 | kappa = 2*sin(phi2-phi_0)/r2; | |
272 | tot_charge = digits[i].charge + digits[j].charge; | |
273 | hist->Fill(kappa,phi_0,tot_charge); | |
b1886074 | 274 | } |
275 | } | |
3ef466c5 | 276 | delete [] digits; |
b1886074 | 277 | } |
278 | ||
4fc9a6a4 | 279 | void AliL3HoughTransformer::TransformLine() |
280 | { | |
281 | //Do a line transform on the data. | |
282 | ||
283 | ||
c52cf5d8 | 284 | AliL3DigitRowData *tempPt = GetDataPointer(); |
285 | if(!tempPt) | |
4fc9a6a4 | 286 | { |
237d3f5c | 287 | LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data") |
288 | <<"No input data "<<ENDLOG; | |
289 | return; | |
290 | } | |
26abc209 | 291 | |
292 | for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++) | |
4fc9a6a4 | 293 | { |
294 | AliL3DigitData *digPt = tempPt->fDigitData; | |
295 | if(i != (Int_t)tempPt->fRow) | |
296 | { | |
297 | printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n"); | |
298 | continue; | |
299 | } | |
300 | for(UInt_t j=0; j<tempPt->fNDigit; j++) | |
301 | { | |
302 | UShort_t charge = digPt[j].fCharge; | |
303 | UChar_t pad = digPt[j].fPad; | |
304 | UShort_t time = digPt[j].fTime; | |
c52cf5d8 | 305 | if(charge < GetThreshold()) |
4fc9a6a4 | 306 | continue; |
307 | Int_t sector,row; | |
308 | Float_t xyz[3]; | |
4ab9f8f0 | 309 | AliL3Transform::Slice2Sector(GetSlice(),i,sector,row); |
310 | AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time); | |
311 | Float_t eta = AliL3Transform::GetEta(xyz); | |
b5a207b4 | 312 | Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice); |
c52cf5d8 | 313 | if(eta_index < 0 || eta_index >= GetNEtaSegments()) |
4fc9a6a4 | 314 | continue; |
315 | ||
316 | //Get the correct histogram: | |
317 | AliL3Histogram *hist = fParamSpace[eta_index]; | |
318 | if(!hist) | |
319 | { | |
320 | printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index); | |
321 | continue; | |
322 | } | |
323 | for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++) | |
324 | { | |
325 | Double_t theta = hist->GetBinCenterX(xbin); | |
326 | Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta); | |
327 | hist->Fill(theta,rho,charge); | |
328 | } | |
329 | } | |
48f3c46f | 330 | AliL3MemHandler::UpdateRowPointer(tempPt); |
4fc9a6a4 | 331 | } |
332 | ||
333 | } |