HLT/hough/AliL3HoughTransformer.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Logging.h"
   9 #include "AliL3HoughTransformer.h"
  10 #include "AliL3MemHandler.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3DigitData.h"
  13 #include "AliL3HistogramAdaptive.h"
  14
  15 #if GCCVERSION == 3
  16 using namespace std;
  17 #endif
  18
  19 //_____________________________________________________________
  20 // AliL3HoughTransformer
  21 //
  22 // Hough transformation class
  23 //
  24
  25 ClassImp(AliL3HoughTransformer)
  26
  27 AliL3HoughTransformer::AliL3HoughTransformer()
  28 {
  29   //Default constructor
  30   fParamSpace = 0;
  31   fDoMC = kFALSE;;
  32 #ifdef do_mc
  33   fTrackID = 0;
  34 #endif
  35 }
  36
  37 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoMC) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
  38 {
  39   //Normal constructor
  40   fParamSpace = 0;
  41   if(DoMC)
  42     {
  43       if(patch==0)
  44         fDoMC = kTRUE;
  45       else
  46         fDoMC = kFALSE;
  47     }
  48 #ifdef do_mc
  49   fTrackID = 0;
  50 #endif
  51 }
  52
  53 AliL3HoughTransformer::~AliL3HoughTransformer()
  54 {
  55   DeleteHistograms();
  56 #ifdef do_mc
  57   if(fTrackID)
  58     {
  59       for(Int_t i=0; i<GetNEtaSegments(); i++)
  60         {
  61           if(!fTrackID[i]) continue;
  62           delete fTrackID[i];
  63         }
  64       delete [] fTrackID;
  65     }
  66 #endif
  67 }
  68
  69 //void AliL3HoughTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
  70
  71 void AliL3HoughTransformer::DeleteHistograms()
  72 {
  73   if(!fParamSpace)
  74     return;
  75   for(Int_t i=0; i<GetNEtaSegments(); i++)
  76     {
  77       if(!fParamSpace[i]) continue;
  78       delete fParamSpace[i];
  79     }
  80   delete [] fParamSpace;
  81 }
  82
  83 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
  84                                              Int_t nybin,Double_t phimin,Double_t phimax)
  85 {
  86   //Create the histograms (parameter space).
  87   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  88   //The arguments give the range and binning;
  89   //nxbin = #bins in kappa
  90   //nybin = #bins in phi0
  91   //pt_min = mimium Pt of track (corresponding to maximum kappa)
  92   //phi_min = mimimum phi0 (degrees)
  93   //phi_max = maximum phi0 (degrees)
  94
  95   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
  96   Double_t torad = AliL3Transform::Pi()/180;
  97   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
  98 }
  99
 100 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
 101                                              Int_t nybin,Double_t ymin,Double_t ymax)
 102 {
 103
 104   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
 105
 106   Char_t histname[256];
 107   for(Int_t i=0; i<GetNEtaSegments(); i++)
 108     {
 109       sprintf(histname,"paramspace_%d",i);
 110       //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.1,2,0.02,nybin,ymin,ymax);
 111       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 112     }
 113
 114 #ifdef do_mc
 115   if(fDoMC)
 116     {
 117       AliL3Histogram *hist = fParamSpace[0];
 118       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 119       cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
 120       fTrackID = new TrackIndex*[GetNEtaSegments()];
 121       for(Int_t i=0; i<GetNEtaSegments(); i++)
 122         fTrackID[i] = new TrackIndex[ncells];
 123     }
 124 #endif
 125 }
 126
 127 void AliL3HoughTransformer::Reset()
 128 {
 129   //Reset all the histograms
 130
 131   if(!fParamSpace)
 132     {
 133       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
 134         <<"No histograms to reset"<<ENDLOG;
 135       return;
 136     }
 137
 138   for(Int_t i=0; i<GetNEtaSegments(); i++)
 139     fParamSpace[i]->Reset();
 140 #ifdef do_mc
 141   if(fDoMC)
 142     {
 143       AliL3Histogram *hist = fParamSpace[0];
 144       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 145       for(Int_t i=0; i<GetNEtaSegments(); i++)
 146         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
 147     }
 148 #endif
 149 }
 150
 151 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
 152 {
 153   //Return the histogram index of the corresponding eta.
 154
 155   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 156   Double_t index = (eta-GetEtaMin())/etaslice;
 157   return (Int_t)index;
 158 }
 159
 160 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
 161 {
 162   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
 163     return 0;
 164   if(!fParamSpace[eta_index])
 165     return 0;
 166   return fParamSpace[eta_index];
 167 }
 168
 169 Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice)
 170 {
 171   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 172   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
 173   return eta;
 174 }
 175
 176 void AliL3HoughTransformer::TransformCircle()
 177 {
 178   //Transform the input data with a circle HT.
 179   //The function loops over all the data, and transforms each pixel with the equations:
 180   //
 181   //kappa = 2/R*sin(phi - phi0)
 182   //
 183   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
 184   //
 185   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
 186   //which histogram in which the pixel should be transformed, the eta-value is calcluated
 187   //and the proper histogram index is found by GetEtaIndex(eta).
 188
 189
 190   AliL3DigitRowData *tempPt = GetDataPointer();
 191   if(!tempPt)
 192     {
 193       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 194         <<"No input data "<<ENDLOG;
 195       return;
 196     }
 197
 198   //Loop over the padrows:
 199   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 200     {
 201       //Get the data on this padrow:
 202       AliL3DigitData *digPt = tempPt->fDigitData;
 203       if(i != (Int_t)tempPt->fRow)
 204         {
 205           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
 206           continue;
 207         }
 208
 209       //Loop over the data on this padrow:
 210       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 211         {
 212           UShort_t charge = digPt[j].fCharge;
 213           UChar_t pad = digPt[j].fPad;
 214           UShort_t time = digPt[j].fTime;
 215           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
 216             continue;
 217
 218           Int_t sector,row;
 219           Float_t xyz[3];
 220
 221           //Transform data to local cartesian coordinates:
 222           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 223           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 224
 225           //Calculate the eta:
 226           Double_t eta = AliL3Transform::GetEta(xyz);
 227
 228           //Get the corresponding index, which determines which histogram to fill:
 229           Int_t eta_index = GetEtaIndex(eta);
 230           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 231             continue;
 232
 233           //Get the correct histogrampointer:
 234           AliL3Histogram *hist = fParamSpace[eta_index];
 235           if(!hist)
 236             {
 237               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
 238               continue;
 239             }
 240
 241           //Do the transformation:
 242           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 243           Float_t phi = AliL3Transform::GetPhi(xyz);
 244
 245           //Fill the histogram along the phirange
 246           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 247             {
 248               Float_t phi0 = hist->GetBinCenterY(b);
 249               Float_t kappa = 2*sin(phi - phi0)/R;
 250               hist->Fill(kappa,phi0,1);//charge);
 251 #ifdef do_mc
 252               if(fDoMC)
 253                 {
 254                   Int_t bin = hist->FindBin(kappa,phi0);
 255                   for(Int_t t=0; t<3; t++)
 256                     {
 257                       Int_t label = digPt[j].fTrackID[t];
 258                       if(label < 0) break;
 259                       UInt_t c;
 260                       for(c=0; c<MaxTrack; c++)
 261                         if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
 262                           break;
 263                       if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 264                       fTrackID[eta_index][bin].fLabel[c] = label;
 265                       fTrackID[eta_index][bin].fNHits[c]++;
 266                     }
 267                 }
 268 #endif
 269             }
 270         }
 271
 272       //Move the data pointer to the next padrow:
 273       AliL3MemHandler::UpdateRowPointer(tempPt);
 274     }
 275 }
 276
 277 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
 278 {
 279   //Circle transform, using combinations of every 2 points lying
 280   //on different padrows and within the same etaslice.
 281
 282   AliL3DigitRowData *tempPt = GetDataPointer();
 283   if(!tempPt)
 284     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 285       <<"No input data "<<ENDLOG;
 286
 287   Int_t counter=0;
 288   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 289     {
 290       counter += tempPt->fNDigit;
 291       AliL3MemHandler::UpdateRowPointer(tempPt);
 292     }
 293
 294   struct Digit {
 295     Int_t row;
 296     Double_t r;
 297     Double_t phi;
 298     Int_t eta_index;
 299     Int_t charge;
 300     Int_t trackID[3];
 301   };
 302
 303   Digit *digits = new Digit[counter];
 304   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
 305
 306   Int_t total_digits=counter;
 307   Int_t sector,row,tot_charge,pad,time,charge;
 308   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
 309   Float_t xyz[3];
 310
 311   counter=0;
 312   tempPt = GetDataPointer();
 313
 314   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 315     {
 316       AliL3DigitData *digPt = tempPt->fDigitData;
 317       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 318         {
 319           charge = digPt[di].fCharge;
 320           pad = digPt[di].fPad;
 321           time = digPt[di].fTime;
 322           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 323           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 324           eta = AliL3Transform::GetEta(xyz);
 325           digits[counter].row = i;
 326           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 327           digits[counter].phi = atan2(xyz[1],xyz[0]);
 328           digits[counter].eta_index = GetEtaIndex(eta);
 329           digits[counter].charge = charge;
 330 #ifdef do_mc
 331           if(fDoMC)
 332             {
 333               digits[counter].trackID[0] = digPt[di].fTrackID[0];
 334               digits[counter].trackID[1] = digPt[di].fTrackID[1];
 335               digits[counter].trackID[2] = digPt[di].fTrackID[2];
 336             }
 337 #endif
 338           counter++;
 339         }
 340       AliL3MemHandler::UpdateRowPointer(tempPt);
 341     }
 342
 343   for(Int_t i=0; i<total_digits; i++)
 344     {
 345       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
 346       Int_t ind = digits[i].eta_index;
 347
 348       for(Int_t j=i+1; j<total_digits; j++)
 349         {
 350           if(digits[i].row == digits[j].row) continue;
 351           if(digits[i].eta_index != digits[j].eta_index) continue;
 352           if(digits[i].row + row_range < digits[j].row) break;
 353
 354           //Get the correct histogrampointer:
 355           AliL3Histogram *hist = fParamSpace[ind];
 356           if(!hist)
 357             {
 358               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
 359               continue;
 360             }
 361
 362           r1 = digits[i].r;
 363           phi1 = digits[i].phi;
 364           r2 = digits[j].r;
 365           phi2 = digits[j].phi;
 366           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
 367           kappa = 2*sin(phi2-phi_0)/r2;
 368           tot_charge = digits[i].charge + digits[j].charge;
 369           hist->Fill(kappa,phi_0,1);//tot_charge);
 370 #ifdef do_mc
 371           if(fDoMC)
 372             {
 373               Int_t bin = hist->FindBin(kappa,phi_0);
 374               for(Int_t l=0; l<3; l++)
 375                 {
 376                   for(Int_t m=0; m<3; m++)
 377                     {
 378                       if(digits[i].trackID[l] == digits[j].trackID[m])
 379                         {
 380                           Int_t label = digits[i].trackID[l];
 381                           if(label < 0) continue;
 382                           UInt_t c;
 383                           for(c=0; c<MaxTrack; c++)
 384                             if(fTrackID[ind][bin].fLabel[c] == label || fTrackID[ind][bin].fNHits[c] == 0)
 385                               break;
 386                           if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
 387                           fTrackID[ind][bin].fLabel[c] = label;
 388                           fTrackID[ind][bin].fNHits[c]++;
 389                         }
 390                     }
 391                 }
 392             }
 393 #endif
 394         }
 395     }
 396   delete [] digits;
 397 }
 398
 399 void AliL3HoughTransformer::TransformLine()
 400 {
 401   //Do a line transform on the data.
 402
 403
 404   AliL3DigitRowData *tempPt = GetDataPointer();
 405   if(!tempPt)
 406     {
 407       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
 408         <<"No input data "<<ENDLOG;
 409       return;
 410     }
 411
 412   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 413     {
 414       AliL3DigitData *digPt = tempPt->fDigitData;
 415       if(i != (Int_t)tempPt->fRow)
 416         {
 417           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
 418           continue;
 419         }
 420       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 421         {
 422           UShort_t charge = digPt[j].fCharge;
 423           UChar_t pad = digPt[j].fPad;
 424           UShort_t time = digPt[j].fTime;
 425           if(charge < GetLowerThreshold())
 426             continue;
 427           Int_t sector,row;
 428           Float_t xyz[3];
 429           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 430           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 431           Float_t eta = AliL3Transform::GetEta(xyz);
 432           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 433           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 434             continue;
 435
 436           //Get the correct histogram:
 437           AliL3Histogram *hist = fParamSpace[eta_index];
 438           if(!hist)
 439             {
 440               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
 441               continue;
 442             }
 443           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
 444             {
 445               Double_t theta = hist->GetBinCenterX(xbin);
 446               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 447               hist->Fill(theta,rho,charge);
 448             }
 449         }
 450       AliL3MemHandler::UpdateRowPointer(tempPt);
 451     }
 452
 453 }
 454
 455 Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
 456 {
 457   if(!fDoMC)
 458     {
 459       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
 460       return -1;
 461     }
 462
 463 #ifdef do_mc
 464   if(eta_index < 0 || eta_index > GetNEtaSegments())
 465     {
 466       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
 467       return -1;
 468     }
 469   AliL3Histogram *hist = fParamSpace[eta_index];
 470   Int_t bin = hist->FindBin(kappa,psi);
 471   Int_t label=-1;
 472   Int_t max=0;
 473   for(UInt_t i=0; i<MaxTrack; i++)
 474     {
 475       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
 476       if(nhits == 0) break;
 477       if(nhits > max)
 478         {
 479           max = nhits;
 480           label = fTrackID[eta_index][bin].fLabel[i];
 481         }
 482     }
 483   return label;
 484 #endif
 485   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
 486   return -1;
 487 }
 488