Made a new abstract base class; AliL3HoughBaseTransformer for different implementations
[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
31
c52cf5d8 32 fParamSpace = 0;
52a2a604 33}
34
4de874d1 35AliL3HoughTransformer::~AliL3HoughTransformer()
36{
99e7186b 37 DeleteHistograms();
4cafa5fc 38}
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
e26acabd 64 Double_t bfact = 0.0029980;
65 Double_t bfield = 0.2;
66 Double_t x = bfact*bfield/pt_min;
b5a207b4 67 CreateHistograms(nxbin,-1.*x,x,nybin,phimin*ToRad,phimax*ToRad);
e26acabd 68}
69
99e7186b 70void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
71 Int_t nybin,Double_t ymin,Double_t ymax)
4cafa5fc 72{
4fc9a6a4 73
c52cf5d8 74 fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
4cafa5fc 75
99e7186b 76 Char_t histname[256];
c52cf5d8 77 for(Int_t i=0; i<GetNEtaSegments(); i++)
4de874d1 78 {
99e7186b 79 sprintf(histname,"paramspace_%d",i);
80 fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
4de874d1 81 }
4cafa5fc 82}
83
e26acabd 84void AliL3HoughTransformer::Reset()
85{
48f3c46f 86 //Reset all the histograms
e26acabd 87
88 if(!fParamSpace)
89 {
90 LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
91 <<"No histograms to reset"<<ENDLOG;
92 return;
93 }
94
c52cf5d8 95 for(Int_t i=0; i<GetNEtaSegments(); i++)
e26acabd 96 fParamSpace[i]->Reset();
e26acabd 97}
98
4de874d1 99
b5a207b4 100Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
101{
3ef466c5 102 //Return the histogram index of the corresponding eta.
103
c52cf5d8 104 Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
105 Double_t index = (eta-GetEtaMin())/etaslice;
b5a207b4 106 return (Int_t)index;
107}
108
4fc9a6a4 109void AliL3HoughTransformer::TransformCircle()
4de874d1 110{
99e7186b 111 //Transform the input data with a circle HT.
3ef466c5 112 //The function loops over all the data, and transforms each pixel with the equations:
113 //
114 //kappa = 2/R*sin(phi - phi0)
115 //
116 //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
117 //
118 //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
119 //which histogram in which the pixel should be transformed, the eta-value is calcluated
120 //and the proper histogram index is found by GetEtaIndex(eta).
121
52a2a604 122
c52cf5d8 123 AliL3DigitRowData *tempPt = GetDataPointer();
124 if(!tempPt)
4de874d1 125 {
4fc9a6a4 126 printf("\nAliL3HoughTransformer::TransformCircle : No input data!!!\n\n");
99e7186b 127 return;
4de874d1 128 }
99e7186b 129
3ef466c5 130 //Loop over the padrows:
c52cf5d8 131 for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; 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:
c52cf5d8 153 fTransform->Slice2Sector(GetSlice(),i,sector,row);
99e7186b 154 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
3ef466c5 155
156 //Calculate the eta:
4fc9a6a4 157 Double_t eta = fTransform->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]);
99e7186b 174 Float_t phi = fTransform->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 190void 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)
197 printf("\nAliL3HoughTransformer::TransformCircleC() : Zero data pointer\n");
198
3ef466c5 199 Int_t counter=0;
c52cf5d8 200 for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
b1886074 201 {
3ef466c5 202 counter += tempPt->fNDigit;
b1886074 203 AliL3MemHandler::UpdateRowPointer(tempPt);
204 }
3ef466c5 205
206 struct Digit {
207 Int_t row;
208 Double_t r;
209 Double_t phi;
210 Int_t eta_index;
211 Int_t charge;
212 };
213
214 Digit *digits = new Digit[counter];
215 cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
216
217 Int_t total_digits=counter;
218 Int_t sector,row,tot_charge,pad,time,charge;
b1886074 219 Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
b1886074 220 Float_t xyz[3];
3ef466c5 221
222 counter=0;
c52cf5d8 223 tempPt = GetDataPointer();
3ef466c5 224
c52cf5d8 225 for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
b1886074 226 {
3ef466c5 227 AliL3DigitData *digPt = tempPt->fDigitData;
228 for(UInt_t di=0; di<tempPt->fNDigit; di++)
b1886074 229 {
3ef466c5 230 charge = digPt[di].fCharge;
b1886074 231 pad = digPt[di].fPad;
232 time = digPt[di].fTime;
c52cf5d8 233 fTransform->Slice2Sector(GetSlice(),i,sector,row);
b1886074 234 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
235 eta = fTransform->GetEta(xyz);
3ef466c5 236 digits[counter].row = i;
237 digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
238 digits[counter].phi = atan2(xyz[1],xyz[0]);
239 digits[counter].eta_index = GetEtaIndex(eta);
240 digits[counter].charge = charge;
241 counter++;
242 }
243 AliL3MemHandler::UpdateRowPointer(tempPt);
244 }
245
246 for(Int_t i=0; i<total_digits; i++)
247 {
c52cf5d8 248 if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
3ef466c5 249 Int_t ind = digits[i].eta_index;
250
251 for(Int_t j=i+1; j<total_digits; j++)
252 {
253 if(digits[i].row == digits[j].row) continue;
254 if(digits[i].eta_index != digits[j].eta_index) continue;
255 if(digits[i].row + row_range < digits[j].row) break;
b1886074 256
257 //Get the correct histogrampointer:
3ef466c5 258 AliL3Histogram *hist = fParamSpace[ind];
b1886074 259 if(!hist)
260 {
3ef466c5 261 printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
b1886074 262 continue;
263 }
3ef466c5 264
265 r1 = digits[i].r;
266 phi1 = digits[i].phi;
267 r2 = digits[j].r;
268 phi2 = digits[j].phi;
269 phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
270 kappa = 2*sin(phi2-phi_0)/r2;
271 tot_charge = digits[i].charge + digits[j].charge;
272 hist->Fill(kappa,phi_0,tot_charge);
b1886074 273 }
274 }
3ef466c5 275 delete [] digits;
b1886074 276}
277
4fc9a6a4 278void AliL3HoughTransformer::TransformLine()
279{
280 //Do a line transform on the data.
281
282
c52cf5d8 283 AliL3DigitRowData *tempPt = GetDataPointer();
284 if(!tempPt)
4fc9a6a4 285 {
286 printf("\nAliL3HoughTransformer::TransformLine : No input data!!!\n\n");
287 return;
288 }
289
c52cf5d8 290 for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
4fc9a6a4 291 {
292 AliL3DigitData *digPt = tempPt->fDigitData;
293 if(i != (Int_t)tempPt->fRow)
294 {
295 printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
296 continue;
297 }
298 for(UInt_t j=0; j<tempPt->fNDigit; j++)
299 {
300 UShort_t charge = digPt[j].fCharge;
301 UChar_t pad = digPt[j].fPad;
302 UShort_t time = digPt[j].fTime;
c52cf5d8 303 if(charge < GetThreshold())
4fc9a6a4 304 continue;
305 Int_t sector,row;
306 Float_t xyz[3];
c52cf5d8 307 fTransform->Slice2Sector(GetSlice(),i,sector,row);
4fc9a6a4 308 fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
309 Float_t eta = fTransform->GetEta(xyz);
b5a207b4 310 Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
c52cf5d8 311 if(eta_index < 0 || eta_index >= GetNEtaSegments())
4fc9a6a4 312 continue;
313
314 //Get the correct histogram:
315 AliL3Histogram *hist = fParamSpace[eta_index];
316 if(!hist)
317 {
318 printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
319 continue;
320 }
321 for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
322 {
323 Double_t theta = hist->GetBinCenterX(xbin);
324 Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
325 hist->Fill(theta,rho,charge);
326 }
327 }
48f3c46f 328 AliL3MemHandler::UpdateRowPointer(tempPt);
4fc9a6a4 329 }
330
331}