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