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