d9dbdd41953693756c0a1b0fba0083764d8e1094
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerGlobal.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliL3StandardIncludes.h"
7
8 #include "AliL3Logging.h"
9 #include "AliL3HoughTransformerGlobal.h"
10 #include "AliL3FileHandler.h"
11 #include "AliL3Transform.h"
12 #include "AliL3DigitData.h"
13 #include "AliL3Track.h"
14 #include "AliL3Histogram.h"
15 #include "AliL3TrackArray.h"
16 #include "AliL3HoughMaxFinder.h"
17 #include "AliL3HoughTrack.h"
18 #include "AliL3SpacePointData.h"
19
20 #include <TH2.h>
21
22 #if GCCVERSION == 3
23 using namespace std;
24 #endif
25
26 //_____________________________________________________________
27 // AliL3HoughTransformerGlobal
28 //
29 // Hough transformation class
30 //
31
32 ClassImp(AliL3HoughTransformerGlobal)
33
34 AliL3HoughTransformerGlobal::AliL3HoughTransformerGlobal()
35 {
36   fPadMin=0;
37   fPadMax=0;
38   fNActiveSlice=0;
39   fMemHandler=0;
40   fEvent=0;
41   fPsi=0;
42   fPtMin=0;
43   fTracks=0;
44   fPeakFinder=0;
45   fSeedPadRow=-1;
46 }
47
48 AliL3HoughTransformerGlobal::AliL3HoughTransformerGlobal(Char_t *path,Int_t event)
49 {
50   strcpy(fPath,path);
51   fEvent=event;
52   fMemHandler=0;
53   fPadMin=0;
54   fPadMax=0;
55   fPsi=0;
56   fPtMin=0;
57   fNActiveSlice=0;
58   fTracks = new AliL3TrackArray("AliL3HoughTrack");
59   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",1000);
60 }
61
62 AliL3HoughTransformerGlobal::~AliL3HoughTransformerGlobal()
63 {
64   if(fPadMin)
65     delete [] fPadMin;
66   if(fPadMax)
67     delete [] fPadMax;
68   if(fTracks)
69     delete fTracks;
70   if(fPeakFinder)
71     delete fPeakFinder;
72   UnloadActiveSlices();
73   
74 }
75
76 void AliL3HoughTransformerGlobal::CreateHistograms(Int_t nxbin,Int_t nybin)
77 {
78   if(fPsi==0)
79     {
80       cerr<<"AliL3HoughTransformerGlobal::CreateHistograms : Call DefineRegion first"<<endl;
81       exit(5);
82     }
83   Float_t psi = fPsi*180./AliL3Transform::Pi();
84   AliL3HoughTransformer::CreateHistograms(nxbin,fPtMin,nybin,-psi,psi);
85 }
86
87 void AliL3HoughTransformerGlobal::TransformCircleC()
88 {
89   
90   Seed *clusters = new Seed[1000];
91   Int_t nclusters = LoadClusterSeeds(clusters);
92   
93   Int_t i=0,sector,row,slice;
94   UInt_t dummy=0;
95   Float_t xyz[3],rpe[3],kappa,psi;
96   
97   for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
98     {
99       if(sl < 0 || sl < 18 && GetSlice() >= 18)
100         slice = sl + 18;
101       else if(sl > 35 || sl > 17 && GetSlice() <= 17)
102         slice = sl-18;
103       else
104         slice = sl;
105       cout<<"Transforming in slice "<<slice<<endl;
106       
107       AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
108       
109       for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
110         {
111           if(padrow == fSeedPadRow) continue;
112           AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
113           
114           for(UInt_t j=0; j<rowPt->fNDigit; j++)
115             {
116               Int_t pad = digits[j].fPad;
117               
118               if(i==1 && pad < fPadMin[padrow])
119                 continue;
120               if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
121                 continue;
122               
123               Int_t time = digits[j].fTime;
124               Int_t charge = digits[j].fCharge;
125               
126               if(charge > GetUpperThreshold() || charge < GetLowerThreshold())
127                 continue;
128               
129               AliL3Transform::Slice2Sector(slice,padrow,sector,row);
130               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
131               
132               Rotate(xyz,i-1-fNActiveSlice);
133               
134               AliL3Transform::XYZtoRPhiEta(rpe,xyz);
135               rpe[0] = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
136               
137               Int_t eta_index = GetEtaIndex(rpe[2]);
138               
139               AliL3Histogram *hist = GetHistogram(eta_index);
140               
141               for(Int_t l=0; l<nclusters; l++)
142                 {
143                   if(clusters[l].fIndex != eta_index) continue;
144                   psi = atan( (clusters[l].fRadius*sin(rpe[1])-rpe[0]*sin(clusters[l].fPhi))/
145                               (clusters[l].fRadius*cos(rpe[1])-rpe[0]*cos(clusters[l].fPhi)) );
146                   kappa = 2*sin(clusters[l].fPhi-psi)/clusters[l].fRadius;
147                   if(fabs(kappa) < 1e-33)
148                     continue;
149                   
150                   hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
151                 }
152             }
153           AliL3MemHandler::UpdateRowPointer(rowPt);
154         }
155     }
156   delete [] clusters;
157 }
158
159 void AliL3HoughTransformerGlobal::TransformCircle()
160 {
161   
162   Int_t i=0,sector,row,slice;
163   UInt_t dummy=0;
164   Float_t xyz[3],rpe[3],kappa,psi;
165   
166   for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
167     {
168       if(sl < 0 || sl < 18 && GetSlice() >= 18)
169         slice = sl + 18;
170       else if(sl > 35 || sl > 17 && GetSlice() <= 17)
171         slice = sl-18;
172       else
173         slice = sl;
174       cout<<"Transforming in slice "<<slice<<endl;
175       
176       AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
177       
178       for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
179         {
180
181           AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
182           
183           for(UInt_t j=0; j<rowPt->fNDigit; j++)
184             {
185               Int_t pad = digits[j].fPad;
186               
187               if(i==1 && pad < fPadMin[padrow])
188                 continue;
189               if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
190                 continue;
191               
192               Int_t time = digits[j].fTime;
193               Int_t charge = digits[j].fCharge;
194               
195               if(charge > GetUpperThreshold() || charge < GetLowerThreshold())
196                 continue;
197               
198               AliL3Transform::Slice2Sector(slice,padrow,sector,row);
199               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
200               
201               Rotate(xyz,i-1-fNActiveSlice);
202               
203               AliL3Transform::XYZtoRPhiEta(rpe,xyz);
204               rpe[0] = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
205               
206               Int_t eta_index = GetEtaIndex(rpe[2]);
207               
208               AliL3Histogram *hist = GetHistogram(eta_index);
209               
210               for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
211                 {
212                   psi = hist->GetBinCenterY(b);
213                   kappa = 2*sin(rpe[1] - psi)/rpe[0];
214                   
215                   if(fabs(kappa) < 1e-33)
216                     continue;
217                   
218                   hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
219                 }
220             }
221           AliL3MemHandler::UpdateRowPointer(rowPt);
222         }
223     }
224 }
225
226 void AliL3HoughTransformerGlobal::VerifyTracks(AliL3TrackArray *tracks,Int_t &index)
227 {
228   for(int i=index; i<tracks->GetNTracks(); i++)
229     {
230       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
231       if(!tr) continue;
232
233       Float_t angle=0;
234       AliL3Transform::Local2GlobalAngle(&angle,GetSlice());
235       if(!tr->CalculateReferencePoint(angle,AliL3Transform::Row2X(fSeedPadRow)))
236         {
237           tracks->Remove(i);
238           continue;
239         }
240       Int_t sector,row;
241       AliL3Transform::Slice2Sector(GetSlice(),fSeedPadRow,sector,row);
242       float xyz[3] = {tr->GetPointX(),tr->GetPointY(),tr->GetPointZ()};
243       AliL3Transform::Global2Raw(xyz,sector,row);
244       if(xyz[1]<0 || xyz[1]>=AliL3Transform::GetNPads(fSeedPadRow) ||
245          xyz[2]<0 || xyz[2]>=AliL3Transform::GetNTimeBins())
246         tracks->Remove(i);
247     }
248   tracks->Compress();
249   index = tracks->GetNTracks();
250 }
251
252 void AliL3HoughTransformerGlobal::FindPeaks(AliL3Histogram *hist,Float_t eta)
253 {
254   
255   fPeakFinder->Reset();
256   fPeakFinder->SetHistogram(hist);
257   fPeakFinder->FindAbsMaxima();
258   if(fPeakFinder->GetWeight(0) < 1000)
259     return;
260   
261   AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks->NextTrack();
262   track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
263   track->SetEta(eta);
264   track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
265 }
266
267 void AliL3HoughTransformerGlobal::Rotate(Float_t *xyz,Int_t rel_slice)
268 {
269   //Rotate coordinates from one slice to the slice relative to it;
270   //-1 means lower
271   //1 means upper
272
273   Float_t angle = (Float_t)rel_slice*AliL3Transform::Deg2Rad(20);
274   Float_t x=xyz[0],y=xyz[1];
275   xyz[0] = x*cos(angle) - y*sin(angle);
276   xyz[1] = x*sin(angle) + y*cos(angle);
277 }
278
279
280 void AliL3HoughTransformerGlobal::DefineRegion(Float_t minpt,Float_t linephi,Int_t seedpadrow)
281 {
282   //Setup the region to be included in the transform
283   //This function should be called only once, since it is the same for all slices.
284   
285   fSeedPadRow=seedpadrow;
286   fPtMin = minpt;
287   
288   fPadMin = new Int_t[AliL3Transform::GetNRows()];
289   fPadMax = new Int_t[AliL3Transform::GetNRows()];
290   
291   //Phirange of data which should be included in the transform
292   //Here we assume a min pt to define it; fPtMin
293   //Based on pt (kappa) and the two points; origo and middle point, we can get psi angle
294   
295   //Calculate the upper point, the lower point is symmetrical around the x-axis
296   //The track in this calculation has a positive charge, and the kappa sign is the opposite:  
297   
298   Float_t xyz[3],phi,phi2;
299   
300   phi = CalculateBorder(xyz,-1); //upward bending track
301   phi2 = CalculateBorder(xyz,1); //downward bending track
302   if(phi2 > phi)
303     phi = phi2;
304   
305   cout<<"Phiangle "<<phi*180/AliL3Transform::Pi()<<" psi "<<fPsi*180/AliL3Transform::Pi()<<" nslices "<<fNActiveSlice<<endl;
306   
307   Float_t rotangle = fNActiveSlice*20;
308   
309   //Calculate the LUT for min/max pad for every padrow, and check which slices we need.
310   Int_t pad,sector,row;
311   for(Int_t i=0; i<AliL3Transform::GetNRows(); i++)
312     {
313       AliL3Transform::Slice2Sector(0,i,sector,row);
314       
315       //Lower boundary:
316       pad = AliL3Transform::GetNPads(i)-1;
317       while(pad >= 0)
318         {
319           AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
320           if(AliL3Transform::GetPhi(xyz) > -1.*phi + AliL3Transform::Deg2Rad(rotangle))
321             fPadMin[i] = pad;
322           else
323             break;
324           pad--;
325         }
326       
327       //Upper boundary
328       pad = 0;
329       while(pad < AliL3Transform::GetNPads(i))
330         {
331           AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
332           if(AliL3Transform::GetPhi(xyz) < phi - AliL3Transform::Deg2Rad(rotangle))
333             fPadMax[i] = pad;
334           else
335             break;
336           pad++;
337         }
338     }
339   
340   cout<<"Padmax "<<fPadMax[155]<<endl;
341
342 }
343
344 Float_t AliL3HoughTransformerGlobal::CalculateBorder(Float_t *xyz,Int_t charge)
345 {
346   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(fSeedPadRow),2) + pow(AliL3Transform::GetMaxY(fSeedPadRow),2));
347
348   Double_t kappa = -charge*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/fPtMin;
349
350   //Calculate the psi angle of the track = emission angle with x-axis in origo
351   if( fabs(lineradius*kappa/2) > 1)
352     {
353       cerr<<"AliL3TransformerGlobal::DefineRegion : Angle too big"<<endl;
354       exit(5);
355     }
356   
357   Int_t nslices=0;
358   Float_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
359   if(charge > 0)
360     fPsi=psi;
361   
362   AliL3Track track;
363   track.SetFirstPoint(0,0,0);
364   track.SetPt(fPtMin);
365   track.SetPsi(psi);
366   track.SetCharge(charge);
367   track.CalculateHelix();
368   
369   Int_t crossingrow=0;
370   if(charge < 0)
371     crossingrow = AliL3Transform::GetNRows()-1;
372   
373   Float_t rotangle;
374   
375  redefine:
376   
377   rotangle = nslices*20;
378   
379   if(nslices==0)//here we are in local slice, so we use the appropriate function. a mess of course...
380     {
381       while(!track.GetCrossingPoint(crossingrow,xyz))
382         {
383           if(crossingrow==0)
384             {     
385               cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point1 on row "<<crossingrow<<endl;
386               exit(5);
387             }
388           crossingrow--;
389         }
390     }
391   else
392     {
393       while(!track.CalculateReferencePoint(AliL3Transform::Deg2Rad(rotangle),AliL3Transform::Row2X(crossingrow)))
394         {
395           if(crossingrow==0)
396             {
397               cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point2 on row "<<crossingrow<<endl;
398               exit(5);
399             }
400           crossingrow--;
401         }
402       
403       xyz[0] = track.GetPointX();
404       xyz[1] = track.GetPointY();
405       xyz[2] = track.GetPointZ();
406     }
407   
408   Float_t phi = atan2(xyz[1],xyz[0]);
409   
410   //Rotate coordinates back to local coordinates:
411   Rotate(xyz,-1*nslices);
412   
413   Int_t sector,row;
414   AliL3Transform::Slice2Sector(0,crossingrow,sector,row);
415   AliL3Transform::Local2Raw(xyz,sector,row);
416   if(xyz[1] < 0)
417     {
418       if(crossingrow>0)
419         {
420           cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Wrong pad, probably a deadzone... "<<xyz[1]<<endl;
421           exit(5);
422         }
423       else
424         return 0;//here you only want the crossing point with the outer padrow
425     }
426   if(xyz[1] >= AliL3Transform::GetNPads(crossingrow)) //Here the range is even one more slice away
427     {
428       cerr<<"One more slice with pad "<<xyz[1]<<endl;
429       nslices++;
430       if(charge < 0)
431         crossingrow = AliL3Transform::GetNRows()-1;
432       goto redefine;
433     }
434   
435   if(nslices > fNActiveSlice)
436     fNActiveSlice = nslices;
437   
438   cout<<"Calculated phi "<<phi<<" and pad "<<xyz[1]<<" on row "<<crossingrow<<endl;
439   return phi;
440 }
441
442 void AliL3HoughTransformerGlobal::LoadActiveSlices(Bool_t binary)
443 {
444   if(fMemHandler)
445     UnloadActiveSlices();
446   fMemHandler = new AliL3MemHandler*[(fNActiveSlice*2 + 1)];
447   Char_t filename[1024];
448   UInt_t dummy;
449   Int_t i=0,slice;
450   for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
451     {
452       if(sl < 0 || sl < 18 && GetSlice() >= 18)
453         slice = sl + 18;
454       else if(sl > 35 || sl > 17 && GetSlice() <= 17)
455         slice = sl-18;
456       else
457         slice = sl;
458       cout<<"Loading slice "<<slice<<endl;
459       fMemHandler[i] = new AliL3FileHandler();
460       fMemHandler[i]->Init(slice,-1);
461       if(binary)
462         {
463           sprintf(filename,"%s/binaries/digits_%d_%d_%d.raw",fPath,fEvent,slice,-1);
464           fMemHandler[i]->SetBinaryInput(filename);
465           fMemHandler[i]->CompBinary2Memory(dummy);
466           fMemHandler[i++]->CloseBinaryInput();
467         }
468       else
469         {
470           fMemHandler[i]->SetAliInput("/data1/AliRoot/1track0.4_digits.root");
471           fMemHandler[i]->AliAltroDigits2Memory(dummy,fEvent);
472           fMemHandler[i++]->CloseAliInput();
473         }
474     }
475   
476 }
477
478 void AliL3HoughTransformerGlobal::UnloadActiveSlices()
479 {
480   if(!fMemHandler)
481     return;
482   for(Int_t i=0; i<=fNActiveSlice*2; i++)
483     {
484       if(!fMemHandler[i]) continue;
485       delete fMemHandler[i];
486     }
487   delete [] fMemHandler;
488   fMemHandler=0;
489 }
490
491 Int_t AliL3HoughTransformerGlobal::LoadClusterSeeds(Seed *seeds)
492 {
493   Char_t filename[1024];
494   UInt_t npoints;
495   sprintf(filename,"%s/hough/points_%d_%d_%d.raw",fPath,fEvent,GetSlice(),-1);
496   //sprintf(filename,"../tracker/points_%d_%d_%d.raw",fEvent,GetSlice(),-1);
497   AliL3MemHandler mem;
498   if(!mem.SetBinaryInput(filename))
499     {
500       cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Cannot open file "<<filename<<endl;
501       exit(5);
502     }
503   AliL3SpacePointData *points = (AliL3SpacePointData*)mem.Allocate();
504   mem.Binary2Memory(npoints,points);
505   mem.CloseBinaryInput();
506   
507   Int_t counter=0;
508   for(UInt_t i=0; i<npoints; i++)
509     {
510       Int_t lrow = (Int_t)points[i].fPadRow;
511       if(lrow < fSeedPadRow) continue;
512       if(lrow > fSeedPadRow) break;
513       if(fabs(points[i].fY) > AliL3Transform::GetMaxY(lrow))
514         {
515           cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Seeds seems not to be in local coordinates: "
516               <<points[i].fY<<endl;
517           exit(5);
518         }
519       Float_t xyz[3] = {AliL3Transform::Row2X(lrow),points[i].fY,points[i].fZ};
520       Double_t eta = AliL3Transform::GetEta(xyz);
521       seeds[counter].fPhi = atan2(xyz[1],xyz[0]);
522       seeds[counter].fRadius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
523       seeds[counter].fIndex = GetEtaIndex(eta);
524       seeds[counter++].fEta = eta;
525     }
526   cout<<"Loaded "<<counter<<" cluster seeds on slice "<<GetSlice()<<" padrow "<<fSeedPadRow<<endl;
527   mem.Free();
528   return counter;
529 }
530
531 void AliL3HoughTransformerGlobal::DisplayActiveRegion(TH2F *hist,Int_t eta_index)
532 {
533   //Fill the active region in a histogram
534   
535   Int_t i=0,sector,row,slice;
536   UInt_t dummy=0;
537   Float_t xyz[3];
538   for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
539     {
540       if(sl < 0 || sl < 18 && GetSlice() >= 18)
541         slice = sl + 18;
542       else if(sl > 35 || sl > 17 && GetSlice() <= 17)
543         slice = sl-18;
544       else
545         slice = sl;
546       cout<<"Displaying slice "<<slice<<endl;
547       AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
548       for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
549         {
550           AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
551           for(UInt_t j=0; j<rowPt->fNDigit; j++)
552             {
553               Int_t pad = digits[j].fPad;
554               if(i==1 && pad < fPadMin[padrow])
555                 continue;
556               if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
557                 continue;
558               Int_t time = digits[j].fTime;
559               AliL3Transform::Slice2Sector(slice,padrow,sector,row);
560               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
561              
562               Rotate(xyz,i-1-fNActiveSlice);
563               
564               Double_t eta = AliL3Transform::GetEta(xyz);
565               if(GetEtaIndex(eta) == eta_index)
566                 hist->Fill(xyz[0],xyz[1]);
567             }
568           AliL3MemHandler::UpdateRowPointer(rowPt);
569         }
570     }
571 }