HLT/hough/AliL3HoughTransformer.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include "AliL3MemHandler.h"
   7 #include "AliL3Logging.h"
   8 #include "AliL3HoughTransformer.h"
   9 #include "AliL3Transform.h"
  10 #include "AliL3DigitData.h"
  11 #include "AliL3Histogram.h"
  12
  13 //_____________________________________________________________
  14 // AliL3HoughTransformer
  15 //
  16 // Hough transformation class
  17 //
  18
  19 ClassImp(AliL3HoughTransformer)
  20
  21 AliL3HoughTransformer::AliL3HoughTransformer()
  22 {
  23   //Default constructor
  24   fParamSpace = 0;
  25 }
  26
  27 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
  28 {
  29   //Normal constructor
  30   fParamSpace = 0;
  31 }
  32
  33 AliL3HoughTransformer::~AliL3HoughTransformer()
  34 {
  35   DeleteHistograms();
  36 }
  37
  38 void AliL3HoughTransformer::DeleteHistograms()
  39 {
  40   if(!fParamSpace)
  41     return;
  42   for(Int_t i=0; i<GetNEtaSegments(); i++)
  43     {
  44       if(!fParamSpace[i]) continue;
  45       delete fParamSpace[i];
  46     }
  47   delete [] fParamSpace;
  48 }
  49
  50 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
  51                                              Int_t nybin,Double_t phimin,Double_t phimax)
  52 {
  53   //Create the histograms (parameter space).
  54   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  55   //The arguments give the range and binning;
  56   //nxbin = #bins in kappa
  57   //nybin = #bins in phi0
  58   //pt_min = mimium Pt of track (corresponding to maximum kappa)
  59   //phi_min = mimimum phi0 (degrees)
  60   //phi_max = maximum phi0 (degrees)
  61
  62   Double_t bfact = 0.0029980;
  63   Double_t bfield = 0.2;
  64   Double_t x = bfact*bfield/pt_min;
  65   Double_t torad = AliL3Transform::Pi()/180;
  66   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
  67 }
  68
  69 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
  70                                              Int_t nybin,Double_t ymin,Double_t ymax)
  71 {
  72
  73   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
  74
  75   Char_t histname[256];
  76   for(Int_t i=0; i<GetNEtaSegments(); i++)
  77     {
  78       sprintf(histname,"paramspace_%d",i);
  79       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
  80     }
  81 }
  82
  83 void AliL3HoughTransformer::Reset()
  84 {
  85   //Reset all the histograms
  86
  87   if(!fParamSpace)
  88     {
  89       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
  90         <<"No histograms to reset"<<ENDLOG;
  91       return;
  92     }
  93
  94   for(Int_t i=0; i<GetNEtaSegments(); i++)
  95     fParamSpace[i]->Reset();
  96 }
  97
  98
  99 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
 100 {
 101   //Return the histogram index of the corresponding eta.
 102
 103   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 104   Double_t index = (eta-GetEtaMin())/etaslice;
 105   return (Int_t)index;
 106 }
 107
 108 void AliL3HoughTransformer::TransformCircle()
 109 {
 110   //Transform the input data with a circle HT.
 111   //The function loops over all the data, and transforms each pixel with the equations:
 112   //
 113   //kappa = 2/R*sin(phi - phi0)
 114   //
 115   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
 116   //
 117   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
 118   //which histogram in which the pixel should be transformed, the eta-value is calcluated
 119   //and the proper histogram index is found by GetEtaIndex(eta).
 120
 121
 122   AliL3DigitRowData *tempPt = GetDataPointer();
 123   if(!tempPt)
 124     {
 125       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 126         <<"No input data "<<ENDLOG;
 127       return;
 128     }
 129
 130   //Loop over the padrows:
 131   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 132     {
 133       //Get the data on this padrow:
 134       AliL3DigitData *digPt = tempPt->fDigitData;
 135       if(i != (Int_t)tempPt->fRow)
 136         {
 137           printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
 138           continue;
 139         }
 140
 141       //Loop over the data on this padrow:
 142       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 143         {
 144           UShort_t charge = digPt[j].fCharge;
 145           UChar_t pad = digPt[j].fPad;
 146           UShort_t time = digPt[j].fTime;
 147           if(charge <= GetThreshold())
 148             continue;
 149           Int_t sector,row;
 150           Float_t xyz[3];
 151
 152           //Transform data to local cartesian coordinates:
 153           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 154           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 155
 156           //Calculate the eta:
 157           Double_t eta = AliL3Transform::GetEta(xyz);
 158
 159           //Get the corresponding index, which determines which histogram to fill:
 160           Int_t eta_index = GetEtaIndex(eta);
 161           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 162             continue;
 163
 164           //Get the correct histogrampointer:
 165           AliL3Histogram *hist = fParamSpace[eta_index];
 166           if(!hist)
 167             {
 168               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
 169               continue;
 170             }
 171
 172           //Do the transformation:
 173           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 174           Float_t phi = AliL3Transform::GetPhi(xyz);
 175
 176           //Fill the histogram along the phirange
 177           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 178             {
 179               Float_t phi0 = hist->GetBinCenterY(b);
 180               Float_t kappa = 2*sin(phi - phi0)/R;
 181               hist->Fill(kappa,phi0,charge);
 182             }
 183         }
 184
 185       //Move the data pointer to the next padrow:
 186       AliL3MemHandler::UpdateRowPointer(tempPt);
 187     }
 188 }
 189
 190 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
 191 {
 192   //Circle transform, using combinations of every 2 points lying
 193   //on different padrows and within the same etaslice.
 194
 195   AliL3DigitRowData *tempPt = GetDataPointer();
 196   if(!tempPt)
 197     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 198       <<"No input data "<<ENDLOG;
 199
 200   Int_t counter=0;
 201   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 202     {
 203       counter += tempPt->fNDigit;
 204       AliL3MemHandler::UpdateRowPointer(tempPt);
 205     }
 206
 207   struct Digit {
 208     Int_t row;
 209     Double_t r;
 210     Double_t phi;
 211     Int_t eta_index;
 212     Int_t charge;
 213   };
 214
 215   Digit *digits = new Digit[counter];
 216   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
 217
 218   Int_t total_digits=counter;
 219   Int_t sector,row,tot_charge,pad,time,charge;
 220   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
 221   Float_t xyz[3];
 222
 223   counter=0;
 224   tempPt = GetDataPointer();
 225
 226   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 227     {
 228       AliL3DigitData *digPt = tempPt->fDigitData;
 229       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 230         {
 231           charge = digPt[di].fCharge;
 232           pad = digPt[di].fPad;
 233           time = digPt[di].fTime;
 234           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 235           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 236           eta = AliL3Transform::GetEta(xyz);
 237           digits[counter].row = i;
 238           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 239           digits[counter].phi = atan2(xyz[1],xyz[0]);
 240           digits[counter].eta_index = GetEtaIndex(eta);
 241           digits[counter].charge = charge;
 242           counter++;
 243         }
 244       AliL3MemHandler::UpdateRowPointer(tempPt);
 245     }
 246
 247   for(Int_t i=0; i<total_digits; i++)
 248     {
 249       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
 250       Int_t ind = digits[i].eta_index;
 251
 252       for(Int_t j=i+1; j<total_digits; j++)
 253         {
 254           if(digits[i].row == digits[j].row) continue;
 255           if(digits[i].eta_index != digits[j].eta_index) continue;
 256           if(digits[i].row + row_range < digits[j].row) break;
 257
 258           //Get the correct histogrampointer:
 259           AliL3Histogram *hist = fParamSpace[ind];
 260           if(!hist)
 261             {
 262               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
 263               continue;
 264             }
 265
 266           r1 = digits[i].r;
 267           phi1 = digits[i].phi;
 268           r2 = digits[j].r;
 269           phi2 = digits[j].phi;
 270           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
 271           kappa = 2*sin(phi2-phi_0)/r2;
 272           tot_charge = digits[i].charge + digits[j].charge;
 273           hist->Fill(kappa,phi_0,tot_charge);
 274         }
 275     }
 276   delete [] digits;
 277 }
 278
 279 void AliL3HoughTransformer::TransformLine()
 280 {
 281   //Do a line transform on the data.
 282
 283
 284   AliL3DigitRowData *tempPt = GetDataPointer();
 285   if(!tempPt)
 286     {
 287       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
 288         <<"No input data "<<ENDLOG;
 289       return;
 290     }
 291
 292   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 293     {
 294       AliL3DigitData *digPt = tempPt->fDigitData;
 295       if(i != (Int_t)tempPt->fRow)
 296         {
 297           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
 298           continue;
 299         }
 300       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 301         {
 302           UShort_t charge = digPt[j].fCharge;
 303           UChar_t pad = digPt[j].fPad;
 304           UShort_t time = digPt[j].fTime;
 305           if(charge < GetThreshold())
 306             continue;
 307           Int_t sector,row;
 308           Float_t xyz[3];
 309           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 310           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 311           Float_t eta = AliL3Transform::GetEta(xyz);
 312           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 313           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 314             continue;
 315
 316           //Get the correct histogram:
 317           AliL3Histogram *hist = fParamSpace[eta_index];
 318           if(!hist)
 319             {
 320               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
 321               continue;
 322             }
 323           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
 324             {
 325               Double_t theta = hist->GetBinCenterX(xbin);
 326               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 327               hist->Fill(theta,rho,charge);
 328             }
 329         }
 330       AliL3MemHandler::UpdateRowPointer(tempPt);
 331     }
 332
 333 }