3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliL3StandardIncludes.h"
8 #include "AliL3Logging.h"
9 #include "AliL3HoughTransformerNew.h"
10 #include "AliL3MemHandler.h"
11 #include "AliL3Transform.h"
12 #include "AliL3DigitData.h"
13 #include "AliL3HistogramAdaptive.h"
14 #include "AliL3TrackArray.h"
15 #include "AliL3HoughTrack.h"
16 #include "AliL3SpacePointData.h"
17 #include "AliL3Vertex.h"
18 #include "AliL3Fitter.h"
24 //_____________________________________________________________
25 // AliL3HoughTransformerNew
27 // Hough transformation class
30 ClassImp(AliL3HoughTransformerNew)
32 AliL3HoughTransformerNew::AliL3HoughTransformerNew()
38 AliL3HoughTransformerNew::AliL3HoughTransformerNew(Int_t slice,Int_t patch,Int_t netasegments) : AliL3HoughTransformer(slice,patch,netasegments)
43 AliL3HoughTransformerNew::~AliL3HoughTransformerNew()
50 void AliL3HoughTransformerNew::CreateHistograms(Int_t nxbins,Float_t xlow,Float_t xup,
51 Int_t nybins,Float_t ylow,Float_t yup,
52 Int_t nzbins,Float_t zlow,Float_t zup)
54 //Create the histogram which contain the hough space
56 sprintf(name,"paramspace_%p",(void*)this);
57 fParamSpace3D = new TH3F(name,"",nxbins,xlow,xup,nybins,ylow,yup,nzbins,zlow,zup);
60 void AliL3HoughTransformerNew::Reset()
63 fParamSpace3D->Reset();
68 void AliL3HoughTransformerNew::TransformLine(Int_t *rowrange,Float_t *phirange)
71 AliL3DigitRowData *tempPt = GetDataPointer();
74 LOG(AliL3Log::kError,"AliL3HoughTransformerNew::TransformLine","Data")
75 <<"No input data "<<ENDLOG;
79 TH3 *hist = fParamSpace3D;
80 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
82 AliL3DigitData *digPt = tempPt->fDigitData;
85 AliL3MemHandler::UpdateRowPointer(tempPt);
88 else if(i > rowrange[1])
90 if(i != (Int_t)tempPt->fRow)
92 printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
95 for(UInt_t j=0; j<tempPt->fNDigit; j++)
97 UShort_t charge = digPt[j].fCharge;
98 UChar_t pad = digPt[j].fPad;
99 UShort_t time = digPt[j].fTime;
100 if(charge < GetLowerThreshold() || charge > GetUpperThreshold()) continue;
104 AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
105 AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
107 Float_t phi = AliL3Transform::GetPhi(xyz);
108 if(phi < phirange[0] || phi > phirange[1])
111 xyz[0] = xyz[0] - AliL3Transform::Row2X(rowrange[0]);
112 Float_t x = xyz[0] + AliL3Transform::Row2X(rowrange[0]);
113 Float_t r = sqrt(x*x + xyz[1]*xyz[1]);
114 Float_t delta = atan(xyz[2]/r);
115 for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
117 Float_t theta = hist->GetXaxis()->GetBinCenter(xbin);
118 Float_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
119 hist->Fill(theta,rho,delta,charge);
122 AliL3MemHandler::UpdateRowPointer(tempPt);
128 Int_t fRow;//padrow index
129 Int_t fCharge;//digits charge
130 Float_t fY;//Y coordinate of the digit
131 Float_t fZ;//Z coordinate of the digit
132 AliL3Digit *fNext;//pointer to the next digit
135 struct AliL3EtaContainer {
136 void *fFirst;//pointer to the first digit in a sequence
137 void *fLast;//pointer to the last digit in a sequence
140 void AliL3HoughTransformerNew::TransformLineC(Int_t *rowrange,Float_t *phirange)
143 AliL3DigitRowData *tempPt = GetDataPointer();
145 LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
146 <<"No input data "<<ENDLOG;
149 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
151 counter += tempPt->fNDigit;
152 AliL3MemHandler::UpdateRowPointer(tempPt);
155 Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
156 AliL3EtaContainer *etaPt = new AliL3EtaContainer[bound];
157 memset(etaPt,0,bound*sizeof(AliL3EtaContainer));
159 AliL3Digit *digits = new AliL3Digit[counter];
160 cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
161 memset(digits,0,counter*sizeof(AliL3Digit));
167 tempPt = GetDataPointer();
169 cout<<"Calculating digits in patch "<<GetPatch()<<endl;
170 for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
172 AliL3DigitData *digPt = tempPt->fDigitData;
173 for(UInt_t di=0; di<tempPt->fNDigit; di++)
175 Int_t charge = digPt[di].fCharge;
176 Int_t pad = digPt[di].fPad;
177 Int_t time = digPt[di].fTime;
178 AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
179 AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
180 Double_t eta = AliL3Transform::GetEta(xyz);
182 Float_t phi = atan2(xyz[1],xyz[0]);
183 if(phi < phirange[0] || phi > phirange[1]) continue;
185 digits[counter].fRow = i;
186 digits[counter].fY = xyz[1];
187 digits[counter].fZ = xyz[2];
188 digits[counter].fCharge = charge;
190 Int_t etaindex = GetEtaIndex(eta);
191 Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
193 if(index > 0 && index < bound)
195 if(etaPt[index].fFirst == 0)
196 etaPt[index].fFirst = (void*)(&digits[counter]);
198 ((AliL3Digit*)(etaPt[index].fLast))->fNext = &digits[counter];
199 etaPt[index].fLast = (void*)(&digits[counter]);
203 AliL3MemHandler::UpdateRowPointer(tempPt);
206 cout<<"Doing the combinatorics"<<endl;
208 AliL3Digit *dPt1,*dPt2;
209 TH3 *hist = fParamSpace3D;
210 for(Int_t e=0; e<GetNEtaSegments(); e++)
212 for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
214 Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
216 for(dPt1 = (AliL3Digit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit*)dPt1->fNext)
218 for(Int_t j=i+1; j<=rowrange[1]; j++)
220 Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
222 for(dPt2 = (AliL3Digit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3Digit*)dPt2->fNext)
224 if(dPt1->fRow == dPt2->fRow)
226 cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
231 Float_t x1 = AliL3Transform::Row2X(dPt1->fRow) - AliL3Transform::Row2X(rowrange[0]);
232 Float_t x2 = AliL3Transform::Row2X(dPt2->fRow) - AliL3Transform::Row2X(rowrange[0]);
233 Float_t y1 = dPt1->fY;
234 Float_t y2 = dPt2->fY;
235 Float_t theta = atan2(x2-x1,y1-y2);
236 Float_t rho = x1*cos(theta)+y1*sin(theta);
237 Float_t r1 = sqrt(pow(AliL3Transform::Row2X(dPt1->fRow),2) + pow(y1,2));
238 Float_t delta1 = atan2(dPt1->fZ,r1);
239 Float_t r2 = sqrt(pow(AliL3Transform::Row2X(dPt2->fRow),2) + pow(y2,2));
240 Float_t delta2 = atan2(dPt2->fZ,r2);
241 Float_t delta = (dPt1->fCharge*delta1+dPt2->fCharge*delta2)/(dPt1->fCharge+dPt2->fCharge);
242 hist->Fill(theta,rho,delta,dPt1->fCharge+dPt2->fCharge);