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