3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
6 #include "AliHLTStandardIncludes.h"
7 #include "AliHLTHoughTest.h"
8 #include "AliHLTModelTrack.h"
9 #include "AliHLTTransform.h"
10 #include "AliHLTHistogram.h"
11 #include "AliHLTTrackArray.h"
12 #include "AliHLTHoughTrack.h"
24 //_____________________________________________________________
28 ClassImp(AliHLTHoughTest)
30 AliHLTHoughTest::AliHLTHoughTest()
37 AliHLTHoughTest::~AliHLTHoughTest()
44 Bool_t AliHLTHoughTest::GenerateTrackData(Double_t pt,Double_t psi,Double_t tgl,Int_t sign,Int_t patch,Int_t minhits)
46 //Generate digits according to given track parameters?
50 fData = new AliHLTSimData[AliHLTTransform::GetNRows(patch)];
51 memset(fData,0,AliHLTTransform::GetNRows(patch)*sizeof(AliHLTSimData));
53 AliHLTModelTrack *track = new AliHLTModelTrack();
58 track->SetCharge(sign);
59 track->SetFirstPoint(0,0,0);
60 track->CalculateHelix();
63 // Int_t temp2[AliHLTTransform::GetNTimeBins()];
64 Int_t * temp2 = new Int_t[AliHLTTransform::GetNTimeBins()];
66 Int_t clustercharge=100;
68 for(Int_t i=AliHLTTransform::GetFirstRow(patch); i<=AliHLTTransform::GetLastRow(patch); i++)
72 if(!track->GetCrossingPoint(i,xyz))
75 Int_t rowindex = i - AliHLTTransform::GetFirstRow(patch);
77 AliHLTTransform::Slice2Sector(0,i,sector,row);
78 AliHLTTransform::Local2Raw(xyz,sector,row);
80 if(xyz[1] < 0 || xyz[1] >= AliHLTTransform::GetNPads(i) || xyz[2] < 0 || xyz[2] >= AliHLTTransform::GetNTimeBins())
83 track->SetPadHit(i,xyz[1]);
84 track->SetTimeHit(i,xyz[2]);
86 Double_t beta = track->GetCrossingAngle(i);
87 track->SetCrossingAngleLUT(i,beta);
88 track->CalculateClusterWidths(i);
90 memset(temp,0,200*sizeof(Int_t));
91 memset(temp2,0,AliHLTTransform::GetNTimeBins()*sizeof(Int_t));
92 Double_t xysigma = sqrt(track->GetParSigmaY2(i));
93 Double_t zsigma = sqrt(track->GetParSigmaZ2(i));
97 for(j=0; j<entries; j++)
99 Int_t pad = TMath::Nint(gRandom->Gaus(xyz[1],xysigma));
100 Int_t time = TMath::Nint(gRandom->Gaus(xyz[2],zsigma));
101 if(pad < 0 || pad >= AliHLTTransform::GetNPads(i) || time < 0 || time >= AliHLTTransform::GetNTimeBins())
113 if(temp[j]==0) continue;
115 Int_t index = j - minpad;
117 if(index < 0 || index >= 10)
119 cerr<<"AliHLTHoughTest::GenerateTrackData : Wrong index "<<index<<endl;
123 Int_t seqcharge = clustercharge*temp[j]/entries;
125 for(Int_t k=0; k<AliHLTTransform::GetNTimeBins(); k++)
127 if(temp2[k]==0) continue;
128 Int_t tindex = k - mintime;
129 if(tindex < 0 || tindex >= 10)
131 cerr<<"AliHLTHoughTest::GenerateTrackData : Wrong timeindex "<<tindex<<" "<<k<<" "<<mintime<<endl;
134 Int_t charge = seqcharge*temp2[k]/entries;
139 //cout<<"row "<<i<<" pad "<<j<<" time "<<k<<" charge "<<charge<<endl;
141 fData[rowindex].fPads[index][tindex]=charge;
144 fData[rowindex].fMinpad=minpad;
145 fData[rowindex].fMintime=mintime;
146 fData[rowindex].fNPads=npads;
150 if(hitcounter < minhits)
155 void AliHLTHoughTest::Transform2Circle(AliHLTHistogram *hist)
160 cerr<<"AliHLTHoughTest::Transform : No data"<<endl;
163 Float_t r,phi,phi0,kappa,xyz[3];
164 Int_t pad,time,charge,sector,row;
165 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
167 Int_t rowindex = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
168 AliHLTTransform::Slice2Sector(0,i,sector,row);
169 for(Int_t j=0; j<fData[rowindex].fNPads; j++)
171 pad = j + fData[rowindex].fMinpad;
172 for(Int_t k=0; k<10; k++)
174 time = k + fData[rowindex].fMintime;
175 charge = fData[rowindex].fPads[j][k];
176 if(charge == 0) continue;
177 AliHLTTransform::Raw2Local(xyz,sector,row,pad,time);
179 r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
181 phi = AliHLTTransform::GetPhi(xyz);
183 for(Int_t k=hist->GetFirstYbin(); k<=hist->GetLastYbin(); k++)
185 phi0 = hist->GetBinCenterY(k);
186 kappa = 2*sin(phi-phi0)/r;
187 hist->Fill(kappa,phi0,charge);
194 void AliHLTHoughTest::Transform2CircleC(AliHLTHistogram *hist)
199 cerr<<"AliHLTHoughTest::TransformC : No data"<<endl;
202 Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
203 Float_t r1,r2,phi1,phi2,phi0,kappa,hit[3],hit2[3];
204 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
206 Int_t rowindex1 = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
207 for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
209 pad1 = d1 + fData[rowindex1].fMinpad;
210 for(Int_t j=0; j<10; j++)
212 time1 = j + fData[rowindex1].fMintime;
213 charge1 = fData[rowindex1].fPads[d1][j];
214 if(charge1==0) continue;
215 AliHLTTransform::Slice2Sector(0,i,sector,row);
216 AliHLTTransform::Raw2Local(hit,sector,row,pad1,time1);
217 r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
218 phi1 = atan2(hit[1],hit[0]);
220 for(Int_t j=i+1; j<=AliHLTTransform::GetLastRow(fCurrentPatch); j++)
222 Int_t rowindex2 = j - AliHLTTransform::GetFirstRow(fCurrentPatch);
223 for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
225 pad2 = d2 + fData[rowindex2].fMinpad;
226 for(Int_t k=0; k<10; k++)
228 time2 = k + fData[rowindex2].fMintime;
229 charge2 = fData[rowindex2].fPads[d2][k];
230 if(charge2==0) continue;
231 AliHLTTransform::Slice2Sector(0,j,sector,row);
232 AliHLTTransform::Raw2Local(hit2,sector,row,pad2,time2);
233 r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
234 phi2 = atan2(hit2[1],hit2[0]);
235 phi0 = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
237 kappa = 2*sin(phi1-phi0) / r1;
238 hist->Fill(kappa,phi0,charge1+charge2);
248 void AliHLTHoughTest::Transform2CircleF(AliHLTHistogram *hist)
250 //Fix one point in the middle of the tpc
254 cerr<<"AliHLTHoughTest::TransformF : No data"<<endl;
257 Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
258 Float_t r1,r2,phi1,phi2,phi0,kappa,hit[3],hit2[3];
260 Int_t rowindex1 = 80 - AliHLTTransform::GetFirstRow(fCurrentPatch);
261 for(Int_t d1=1; d1<fData[rowindex1].fNPads; d1++)
263 pad1 = d1 + fData[rowindex1].fMinpad;
265 AliHLTTransform::Slice2Sector(0,80,sector,row);
266 AliHLTTransform::Raw2Local(hit,sector,row,pad1,0);
267 r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
268 phi1 = atan2(hit[1],hit[0]);
270 for(Int_t j=0; j<10; j++)
272 time1 = j + fData[rowindex1].fMintime;
273 charge1 = fData[rowindex1].fPads[d1][j];
274 if(charge1==0) continue;
276 for(Int_t j=0; j<=AliHLTTransform::GetLastRow(fCurrentPatch); j++)
279 Int_t rowindex2 = j - AliHLTTransform::GetFirstRow(fCurrentPatch);
280 for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
282 pad2 = d2 + fData[rowindex2].fMinpad;
283 for(Int_t k=0; k<10; k++)
285 time2 = k + fData[rowindex2].fMintime;
286 charge2 = fData[rowindex2].fPads[d2][k];
287 if(charge2==0) continue;
288 AliHLTTransform::Slice2Sector(0,j,sector,row);
289 AliHLTTransform::Raw2Local(hit2,sector,row,pad2,time2);
290 r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
291 phi2 = atan2(hit2[1],hit2[0]);
292 phi0 = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
294 kappa = 2*sin(phi1-phi0) / r1;
295 hist->Fill(kappa,phi0,charge1+charge2);
296 //cout<<"Filling "<<kappa<<" psi "<<phi0<<" charge "<<charge1<<endl;
305 void AliHLTHoughTest::Transform2Line(AliHLTHistogram *hist,Int_t *rowrange)
310 cerr<<"AliHLTHoughTest::Transform2Line : No data"<<endl;
314 Int_t pad,time,charge,sector,row;
315 Float_t hit[3],theta,rho;
316 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
322 Int_t rowindex = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
323 for(Int_t d=0; d<fData[rowindex].fNPads; d++)
325 pad = d + fData[rowindex].fMinpad;
326 for(Int_t j=0; j<10; j++)
328 time = j + fData[rowindex].fMintime;
329 charge = fData[rowindex].fPads[d][j];
330 if(charge==0) continue;
331 AliHLTTransform::Slice2Sector(0,i,sector,row);
332 AliHLTTransform::Raw2Local(hit,sector,row,pad,time);
334 hit[0] = hit[0] - AliHLTTransform::Row2X(rowrange[0]);
336 for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
338 theta = hist->GetBinCenterX(xbin);
339 rho = hit[0]*cos(theta) + hit[1]*sin(theta);
340 hist->Fill(theta,rho,charge);
347 void AliHLTHoughTest::Transform2LineC(AliHLTHistogram *hist,Int_t *rowrange)
352 cerr<<"AliHLTHoughTest::Transform2Line : No data"<<endl;
356 Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
357 Float_t theta,rho,hit[3],hit2[3];
358 for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
360 Int_t rowindex1 = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
361 for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
363 pad1 = d1 + fData[rowindex1].fMinpad;
364 for(Int_t j=0; j<10; j++)
366 time1 = j + fData[rowindex1].fMintime;
367 charge1 = fData[rowindex1].fPads[d1][j];
368 if(charge1==0) continue;
369 AliHLTTransform::Slice2Sector(0,i,sector,row);
370 AliHLTTransform::Raw2Local(hit,sector,row,pad1,time1);
372 hit[0] = hit[0] - AliHLTTransform::Row2X(rowrange[0]);
374 for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
376 Int_t rowindex2 = i2 - AliHLTTransform::GetFirstRow(fCurrentPatch);
377 for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
379 pad2 = d2 + fData[rowindex2].fMinpad;
380 for(Int_t k=0; k<10; k++)
382 time2 = k + fData[rowindex2].fMintime;
383 charge2 = fData[rowindex2].fPads[d2][k];
384 if(charge2==0) continue;
385 AliHLTTransform::Slice2Sector(0,i2,sector,row);
386 AliHLTTransform::Raw2Local(hit2,sector,row,pad2,time2);
388 hit2[0] = hit2[0] - AliHLTTransform::Row2X(rowrange[0]);
390 theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
391 rho = hit[0]*cos(theta)+hit[1]*sin(theta);
392 hist->Fill(theta,rho,charge1+charge2);
401 void AliHLTHoughTest::FillImage(TH2 *hist,Int_t displayrow)
406 cerr<<"AliHLTHoughTest::FillImage : No data to fill"<<endl;
410 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
412 Int_t rowindex = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
414 if(i != displayrow) continue;
416 //cout<<"row "<<i<<" fNPads "<<fData[rowindex].fNPads<<endl;
417 for(Int_t j=0; j<fData[rowindex].fNPads; j++)
419 Int_t pad = j + fData[rowindex].fMinpad;
420 for(Int_t k=0; k<10; k++)
422 Int_t time = k + fData[rowindex].fMintime;
423 Int_t charge = fData[rowindex].fPads[j][k];
424 if(charge==0) continue;
425 //cout<<i<<" "<<pad<<" "<<time<<" "<<charge<<endl;
428 AliHLTTransform::Slice2Sector(0,i,sector,row);
429 AliHLTTransform::Raw2Local(xyz,sector,row,pad,time);
431 hist->Fill(pad,time,charge);
433 hist->Fill(xyz[0],xyz[1],charge);
441 void AliHLTHoughTest::Transform2Line3D(TH3 *hist,Int_t *rowrange,Float_t *phirange)
446 cerr<<"AliHLTHoughTest::Transform2Line : No data"<<endl;
450 Int_t pad,time,charge,sector,row;
451 Float_t hit[3],theta,rho,r,delta;
452 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
460 Int_t rowindex = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
461 for(Int_t d=0; d<fData[rowindex].fNPads; d++)
463 pad = d + fData[rowindex].fMinpad;
464 for(Int_t j=0; j<10; j++)
466 time = j + fData[rowindex].fMintime;
467 charge = fData[rowindex].fPads[d][j];
468 if(charge==0) continue;
469 AliHLTTransform::Slice2Sector(0,i,sector,row);
470 AliHLTTransform::Raw2Local(hit,sector,row,pad,time);
472 Float_t phi = AliHLTTransform::GetPhi(hit);
473 if(phi < phirange[0] || phi > phirange[1])
476 hit[0] = hit[0] - AliHLTTransform::Row2X(rowrange[0]);
477 Float_t x = hit[0] + AliHLTTransform::Row2X(rowrange[0]);
478 r = sqrt(x*x + hit[1]*hit[1]);
479 delta = atan(hit[2]/r);
481 for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
483 theta = hist->GetXaxis()->GetBinCenter(xbin);
484 rho = hit[0]*cos(theta) + hit[1]*sin(theta);
485 hist->Fill(theta,rho,delta,charge);
492 void AliHLTHoughTest::Transform2LineC3D(TH3 *hist,Int_t *rowrange)
497 cerr<<"AliHLTHoughTest::Transform2Line : No data"<<endl;
501 Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
502 Float_t theta,rho,hit[3],hit2[3],r1,r2,delta,delta1,delta2;
503 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
509 Int_t rowindex1 = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
510 for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
512 pad1 = d1 + fData[rowindex1].fMinpad;
513 for(Int_t j=0; j<10; j++)
515 time1 = j + fData[rowindex1].fMintime;
516 charge1 = fData[rowindex1].fPads[d1][j];
517 if(charge1==0) continue;
518 AliHLTTransform::Slice2Sector(0,i,sector,row);
519 AliHLTTransform::Raw2Local(hit,sector,row,pad1,time1);
520 r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
521 delta1 = atan(hit[2]/r1);
522 hit[0] = hit[0] - AliHLTTransform::Row2X(rowrange[0]);
524 for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
526 Int_t rowindex2 = i2 - AliHLTTransform::GetFirstRow(fCurrentPatch);
527 for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
529 pad2 = d2 + fData[rowindex2].fMinpad;
530 for(Int_t k=0; k<10; k++)
532 time2 = k + fData[rowindex2].fMintime;
533 charge2 = fData[rowindex2].fPads[d2][k];
534 if(charge2==0) continue;
535 AliHLTTransform::Slice2Sector(0,i2,sector,row);
536 AliHLTTransform::Raw2Local(hit2,sector,row,pad2,time2);
537 r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
538 delta2 = atan(hit2[2]/r2);
539 delta = (charge1*delta1 + charge2*delta2)/(charge1+charge2);
540 hit2[0] = hit2[0] - AliHLTTransform::Row2X(rowrange[0]);
542 theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
543 rho = hit[0]*cos(theta)+hit[1]*sin(theta);
544 hist->Fill(theta,rho,delta,charge1+charge2);
553 void AliHLTHoughTest::TransformLines2Circle(TH3 *hist,AliHLTTrackArray *tracks)
556 for(Int_t i=0; i<tracks->GetNTracks(); i++)
558 AliHLTHoughTrack *tr = (AliHLTHoughTrack*)tracks->GetCheckedTrack(i);
560 Int_t middlerow = (Int_t)(tr->GetFirstRow()+(tr->GetLastRow()-tr->GetFirstRow())/2);
562 tr->GetLineCrossingPoint(middlerow,hit);
563 hit[0] += AliHLTTransform::Row2X(tr->GetFirstRow());
564 Float_t r = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
565 hit[2] = r*tr->GetTgl();
566 Float_t phi = atan2(hit[1],hit[0]);
567 Float_t theta = tr->GetPsiLine() - AliHLTTransform::Pi()/2;
568 Float_t psi = 2*phi - theta;
569 Float_t kappa = 2/r*sin(phi-psi);
570 Float_t delta = atan(hit[2]/r);
571 hist->Fill(kappa,psi,delta,tr->GetWeight());
576 void AliHLTHoughTest::Transform2Center(AliHLTHistogram *hist)
578 //Choose parameter space to be center of curvature.
582 cerr<<"AliHLTHoughTest::TransformC : No data"<<endl;
585 Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
586 Float_t r1,r2,phi1,phi2,hit[3],hit2[3];
587 //Float_t phi_0,kappa;
588 for(Int_t i=AliHLTTransform::GetFirstRow(fCurrentPatch); i<=AliHLTTransform::GetLastRow(fCurrentPatch); i++)
590 Int_t rowindex1 = i - AliHLTTransform::GetFirstRow(fCurrentPatch);
591 for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
593 pad1 = d1 + fData[rowindex1].fMinpad;
594 for(Int_t j=0; j<10; j++)
596 time1 = j + fData[rowindex1].fMintime;
597 charge1 = fData[rowindex1].fPads[d1][j];
598 if(charge1==0) continue;
599 AliHLTTransform::Slice2Sector(0,i,sector,row);
600 AliHLTTransform::Raw2Local(hit,sector,row,pad1,time1);
601 r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1])/2;
602 phi1 = atan2(hit[1],hit[0]);
604 for(Int_t j=i+1; j<=AliHLTTransform::GetLastRow(fCurrentPatch); j++)
606 Int_t rowindex2 = j - AliHLTTransform::GetFirstRow(fCurrentPatch);
607 for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
609 pad2 = d2 + fData[rowindex2].fMinpad;
610 for(Int_t k=0; k<10; k++)
612 time2 = k + fData[rowindex2].fMintime;
613 charge2 = fData[rowindex2].fPads[d2][k];
614 if(charge2==0) continue;
615 AliHLTTransform::Slice2Sector(0,j,sector,row);
616 AliHLTTransform::Raw2Local(hit2,sector,row,pad2,time2);
617 r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1])/2;
618 phi2 = atan2(hit2[1],hit2[0]);
619 Float_t yc = (r2-(r1/cos(phi1))*cos(phi2))/(sin(phi2)-tan(phi1)*cos(phi2));
620 Float_t xc = r1/cos(phi1) - (r2*tan(phi1)-r1*sin(phi1)*cos(phi2))/(sin(phi2)-tan(phi1)*cos(phi2));
621 hist->Fill(xc,yc,charge1+charge2);
631 void AliHLTHoughTest::FindAbsMaxima(TH3 *hist,Int_t zsearch,Float_t &maxx,Float_t &maxy,Float_t &maxz,Int_t &maxvalue) const
633 //Find peaks in the Hough space
634 TH1 *h1 = hist->Project3D("z");
636 TAxis *z = hist->GetZaxis();
638 Int_t n=0,i,zbin[50];
639 for(i=z->GetFirst()+1; i<z->GetLast()-1; i++)
641 int bin1 = h1->GetBin(i-1);
642 int bin2 = h1->GetBin(i);
643 int bin3 = h1->GetBin(i+1);
644 if(h1->GetBinContent(bin2) > h1->GetBinContent(bin1) && h1->GetBinContent(bin2) > h1->GetBinContent(bin3))
647 zpeak[n++]=h1->GetBinCenter(bin2);
651 Int_t zrange[2] = {z->GetFirst(),z->GetLast()};
652 z->SetRange(zbin[0]-zsearch,zbin[0]+zsearch);
654 TH1 *h2 = hist->Project3D("yx");
655 z->SetRange(zrange[0],zrange[1]);
657 TAxis *x = h2->GetXaxis();
658 TAxis *y = h2->GetYaxis();
662 for(Int_t xbin=x->GetFirst(); xbin<=x->GetLast(); xbin++)
664 Float_t xvalue = x->GetBinCenter(xbin);
665 for(Int_t ybin=y->GetFirst(); ybin<=y->GetLast(); ybin++)
667 Float_t yvalue = y->GetBinCenter(ybin);
668 Int_t value = (Int_t)h2->GetBinContent(xbin,ybin);