Some additional changes related to the previous changes. AliL3Transform
[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 "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 20ClassImp(AliL3HoughTransformer)
21
22AliL3HoughTransformer::AliL3HoughTransformer()
23{
24 //Default constructor
c52cf5d8 25 fParamSpace = 0;
4de874d1 26}
27
c52cf5d8 28AliL3HoughTransformer::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 34AliL3HoughTransformer::~AliL3HoughTransformer()
35{
99e7186b 36 DeleteHistograms();
4cafa5fc 37}
38
99e7186b 39void 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 51void 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 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
4de874d1 98
b5a207b4 99Int_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 108void 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 196void 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 292void 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}