Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[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(Float_t ptmin,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   AliL3HoughTransformer::CreateHistograms(nxbin,fPtMin,nybin,-fPsi,fPsi);
84   //AliL3HoughTransformer::CreateHistograms(fPtMin,ptmax,ptres,nybin,fPsi);
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                   //hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
150                   hist->Fill(kappa,psi,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,charge);
219                   //hist->Fill(kappa,psi,(int)rint(log((Float_t)charge)));
220                 }
221             }
222           AliL3MemHandler::UpdateRowPointer(rowPt);
223         }
224     }
225 }
226
227 void AliL3HoughTransformerGlobal::VerifyTracks(AliL3TrackArray *tracks,Int_t &index)
228 {
229   for(int i=index; i<tracks->GetNTracks(); i++)
230     {
231       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
232       if(!tr) continue;
233
234       Float_t angle=0;
235       AliL3Transform::Local2GlobalAngle(&angle,GetSlice());
236       if(!tr->CalculateReferencePoint(angle,AliL3Transform::Row2X(fSeedPadRow)))
237         {
238           tracks->Remove(i);
239           continue;
240         }
241       Int_t sector,row;
242       AliL3Transform::Slice2Sector(GetSlice(),fSeedPadRow,sector,row);
243       float xyz[3] = {tr->GetPointX(),tr->GetPointY(),tr->GetPointZ()};
244       AliL3Transform::Global2Raw(xyz,sector,row);
245       if(xyz[1]<0 || xyz[1]>=AliL3Transform::GetNPads(fSeedPadRow) ||
246          xyz[2]<0 || xyz[2]>=AliL3Transform::GetNTimeBins())
247         tracks->Remove(i);
248     }
249   tracks->Compress();
250   index = tracks->GetNTracks();
251 }
252
253 void AliL3HoughTransformerGlobal::FindPeaks(AliL3Histogram *hist,Float_t eta)
254 {
255   
256   fPeakFinder->Reset();
257   fPeakFinder->SetHistogram(hist);
258   fPeakFinder->FindAbsMaxima();
259   if(fPeakFinder->GetWeight(0) < 1000)
260     return;
261   
262   AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks->NextTrack();
263   track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
264   track->SetEta(eta);
265   track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
266 }
267
268 void AliL3HoughTransformerGlobal::Rotate(Float_t *xyz,Int_t rel_slice)
269 {
270   //Rotate coordinates from one slice to the slice relative to it;
271   //-1 means lower
272   //1 means upper
273
274   Float_t angle = (Float_t)rel_slice*AliL3Transform::Deg2Rad(20);
275   Float_t x=xyz[0],y=xyz[1];
276   xyz[0] = x*cos(angle) - y*sin(angle);
277   xyz[1] = x*sin(angle) + y*cos(angle);
278 }
279
280
281 void AliL3HoughTransformerGlobal::DefineRegion(Float_t minpt,Float_t linephi,Int_t seedpadrow)
282 {
283   //Setup the region to be included in the transform
284   //This function should be called only once, since it is the same for all slices.
285   
286   fSeedPadRow=seedpadrow;
287   fPtMin = minpt;
288   
289   fPadMin = new Int_t[AliL3Transform::GetNRows()];
290   fPadMax = new Int_t[AliL3Transform::GetNRows()];
291   
292   //Phirange of data which should be included in the transform
293   //Here we assume a min pt to define it; fPtMin
294   //Based on pt (kappa) and the two points; origo and middle point, we can get psi angle
295   
296   //Calculate the upper point, the lower point is symmetrical around the x-axis
297   //The track in this calculation has a positive charge, and the kappa sign is the opposite:  
298   
299   Float_t xyz[3],phi,phi2;
300   
301   phi = CalculateBorder(xyz,-1); //upward bending track
302   phi2 = CalculateBorder(xyz,1); //downward bending track
303   if(phi2 > phi)
304     phi = phi2;
305   
306   cout<<"Phiangle "<<phi*180/AliL3Transform::Pi()<<" psi "<<fPsi*180/AliL3Transform::Pi()<<" nslices "<<fNActiveSlice<<endl;
307   
308   Float_t rotangle = fNActiveSlice*20;
309   
310   //Calculate the LUT for min/max pad for every padrow, and check which slices we need.
311   Int_t pad,sector,row;
312   for(Int_t i=0; i<AliL3Transform::GetNRows(); i++)
313     {
314       AliL3Transform::Slice2Sector(0,i,sector,row);
315       
316       //Lower boundary:
317       pad = AliL3Transform::GetNPads(i)-1;
318       while(pad >= 0)
319         {
320           AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
321           if(AliL3Transform::GetPhi(xyz) > -1.*phi + AliL3Transform::Deg2Rad(rotangle))
322             fPadMin[i] = pad;
323           else
324             break;
325           pad--;
326         }
327       
328       //Upper boundary
329       pad = 0;
330       while(pad < AliL3Transform::GetNPads(i))
331         {
332           AliL3Transform::Raw2Local(xyz,sector,row,pad,0);
333           if(AliL3Transform::GetPhi(xyz) < phi - AliL3Transform::Deg2Rad(rotangle))
334             fPadMax[i] = pad;
335           else
336             break;
337           pad++;
338         }
339     }
340   
341   cout<<"Padmax "<<fPadMax[155]<<endl;
342
343 }
344
345 Float_t AliL3HoughTransformerGlobal::CalculateBorder(Float_t *xyz,Int_t charge)
346 {
347   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(fSeedPadRow),2) + pow(AliL3Transform::GetMaxY(fSeedPadRow),2));
348
349   Double_t kappa = -charge*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/fPtMin;
350
351   //Calculate the psi angle of the track = emission angle with x-axis in origo
352   if( fabs(lineradius*kappa/2) > 1)
353     {
354       cerr<<"AliL3TransformerGlobal::DefineRegion : Angle too big"<<endl;
355       exit(5);
356     }
357   
358   Int_t nslices=0;
359   Float_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
360   if(charge > 0)
361     fPsi=psi;
362   
363   cout<<"Calculated psi-angle "<<fPsi<<endl;
364   
365   AliL3Track track;
366   track.SetFirstPoint(0,0,0);
367   track.SetPt(fPtMin);
368   track.SetPsi(psi);
369   track.SetCharge(charge);
370   track.CalculateHelix();
371   
372   Int_t crossingrow=0;
373   if(charge < 0)
374     crossingrow = AliL3Transform::GetNRows()-1;
375   
376   Float_t rotangle;
377   
378  redefine:
379   
380   rotangle = nslices*20;
381   
382   if(nslices==0)//here we are in local slice, so we use the appropriate function. a mess of course...
383     {
384       while(!track.GetCrossingPoint(crossingrow,xyz))
385         {
386           if(crossingrow==0)
387             {     
388               cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point1 on row "<<crossingrow<<endl;
389               exit(5);
390             }
391           crossingrow--;
392         }
393     }
394   else
395     {
396       while(!track.CalculateReferencePoint(AliL3Transform::Deg2Rad(rotangle),AliL3Transform::Row2X(crossingrow)))
397         {
398           if(crossingrow==0)
399             {
400               cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Error calculating point2 on row "<<crossingrow<<endl;
401               exit(5);
402             }
403           crossingrow--;
404         }
405       
406       xyz[0] = track.GetPointX();
407       xyz[1] = track.GetPointY();
408       xyz[2] = track.GetPointZ();
409     }
410   
411   Float_t phi = atan2(xyz[1],xyz[0]);
412   
413   //Rotate coordinates back to local coordinates:
414   Rotate(xyz,-1*nslices);
415   
416   Int_t sector,row;
417   AliL3Transform::Slice2Sector(0,crossingrow,sector,row);
418   AliL3Transform::Local2Raw(xyz,sector,row);
419   if(xyz[1] < 0)
420     {
421       if(crossingrow>0)
422         {
423           cerr<<"AliL3HoughTransformerGlobal::DefineRegion : Wrong pad, probably a deadzone... "<<xyz[1]<<endl;
424           exit(5);
425         }
426       else
427         return 0;//here you only want the crossing point with the outer padrow
428     }
429   if(xyz[1] >= AliL3Transform::GetNPads(crossingrow)) //Here the range is even one more slice away
430     {
431       cerr<<"One more slice with pad "<<xyz[1]<<endl;
432       nslices++;
433       if(charge < 0)
434         crossingrow = AliL3Transform::GetNRows()-1;
435       goto redefine;
436     }
437   
438   if(nslices > fNActiveSlice)
439     fNActiveSlice = nslices;
440   
441   cout<<"Calculated phi "<<phi<<" and pad "<<xyz[1]<<" on row "<<crossingrow<<endl;
442   return phi;
443 }
444
445 void AliL3HoughTransformerGlobal::LoadActiveSlices(Bool_t binary)
446 {
447   if(fMemHandler)
448     UnloadActiveSlices();
449   fMemHandler = new AliL3MemHandler*[(fNActiveSlice*2 + 1)];
450   Char_t filename[1024];
451   UInt_t dummy;
452   Int_t i=0,slice;
453   for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
454     {
455       if(sl < 0 || sl < 18 && GetSlice() >= 18)
456         slice = sl + 18;
457       else if(sl > 35 || sl > 17 && GetSlice() <= 17)
458         slice = sl-18;
459       else
460         slice = sl;
461       cout<<"Loading slice "<<slice<<endl;
462       fMemHandler[i] = new AliL3FileHandler();
463       fMemHandler[i]->Init(slice,-1);
464       if(binary)
465         {
466           sprintf(filename,"%s/binaries/digits_%d_%d_%d.raw",fPath,fEvent,slice,-1);
467           fMemHandler[i]->SetBinaryInput(filename);
468           fMemHandler[i]->CompBinary2Memory(dummy);
469           fMemHandler[i++]->CloseBinaryInput();
470         }
471       else
472         {
473           sprintf(filename,"%s/digitfile.root",fPath);
474           fMemHandler[i]->SetAliInput(filename);
475           fMemHandler[i]->AliAltroDigits2Memory(dummy,fEvent);
476           fMemHandler[i++]->CloseAliInput();
477         }
478     }
479   
480 }
481
482 void AliL3HoughTransformerGlobal::UnloadActiveSlices()
483 {
484   if(!fMemHandler)
485     return;
486   for(Int_t i=0; i<=fNActiveSlice*2; i++)
487     {
488       if(!fMemHandler[i]) continue;
489       delete fMemHandler[i];
490     }
491   delete [] fMemHandler;
492   fMemHandler=0;
493 }
494
495 Int_t AliL3HoughTransformerGlobal::LoadClusterSeeds(Seed *seeds)
496 {
497   Char_t filename[1024];
498   UInt_t npoints;
499   sprintf(filename,"%s/hough/points_%d_%d_%d.raw",fPath,fEvent,GetSlice(),-1);
500   //sprintf(filename,"../tracker/points_%d_%d_%d.raw",fEvent,GetSlice(),-1);
501   AliL3MemHandler mem;
502   if(!mem.SetBinaryInput(filename))
503     {
504       cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Cannot open file "<<filename<<endl;
505       exit(5);
506     }
507   AliL3SpacePointData *points = (AliL3SpacePointData*)mem.Allocate();
508   mem.Binary2Memory(npoints,points);
509   mem.CloseBinaryInput();
510   
511   Int_t counter=0;
512   for(UInt_t i=0; i<npoints; i++)
513     {
514       Int_t lrow = (Int_t)points[i].fPadRow;
515       if(lrow < fSeedPadRow) continue;
516       if(lrow > fSeedPadRow) break;
517       if(fabs(points[i].fY) > AliL3Transform::GetMaxY(lrow))
518         {
519           cerr<<"AliL3HoughTransformerGlobal::LoadClusterSeeds : Seeds seems not to be in local coordinates: "
520               <<points[i].fY<<endl;
521           exit(5);
522         }
523       Float_t xyz[3] = {AliL3Transform::Row2X(lrow),points[i].fY,points[i].fZ};
524       Double_t eta = AliL3Transform::GetEta(xyz);
525       seeds[counter].fPhi = atan2(xyz[1],xyz[0]);
526       seeds[counter].fRadius = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
527       seeds[counter].fIndex = GetEtaIndex(eta);
528       seeds[counter++].fEta = eta;
529     }
530   cout<<"Loaded "<<counter<<" cluster seeds on slice "<<GetSlice()<<" padrow "<<fSeedPadRow<<endl;
531   mem.Free();
532   return counter;
533 }
534
535 void AliL3HoughTransformerGlobal::DisplayActiveRegion(TH2F *hist,Int_t eta_index)
536 {
537   //Fill the active region in a histogram
538   
539   Int_t i=0,sector,row,slice;
540   UInt_t dummy=0;
541   Float_t xyz[3];
542   for(Int_t sl=GetSlice()-fNActiveSlice; sl<=GetSlice()+fNActiveSlice; sl++)
543     {
544       if(sl < 0 || sl < 18 && GetSlice() >= 18)
545         slice = sl + 18;
546       else if(sl > 35 || sl > 17 && GetSlice() <= 17)
547         slice = sl-18;
548       else
549         slice = sl;
550       cout<<"Displaying slice "<<slice<<endl;
551       AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fMemHandler[i++]->GetDataPointer(dummy);
552       for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
553         {
554           AliL3DigitData *digits = (AliL3DigitData*)rowPt->fDigitData;
555           for(UInt_t j=0; j<rowPt->fNDigit; j++)
556             {
557               Int_t pad = digits[j].fPad;
558               if(i==1 && pad < fPadMin[padrow])
559                 continue;
560               if(i==fNActiveSlice*2+1 && pad>fPadMax[padrow])
561                 continue;
562               Int_t time = digits[j].fTime;
563               AliL3Transform::Slice2Sector(slice,padrow,sector,row);
564               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
565              
566               Rotate(xyz,i-1-fNActiveSlice);
567               
568               Double_t eta = AliL3Transform::GetEta(xyz);
569               if(GetEtaIndex(eta) == eta_index)
570                 hist->Fill(xyz[0],xyz[1]);
571             }
572           AliL3MemHandler::UpdateRowPointer(rowPt);
573         }
574     }
575 }