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 }