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   fParamSpace = 0;
  26 }
  27
  28 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
  29 {
  30   //Normal constructor
  31   fParamSpace = 0;
  32 }
  33
  34 AliL3HoughTransformer::~AliL3HoughTransformer()
  35 {
  36   DeleteHistograms();
  37 }
  38
  39 void AliL3HoughTransformer::DeleteHistograms()
  40 {
  41   if(!fParamSpace)
  42     return;
  43   for(Int_t i=0; i<GetNEtaSegments(); i++)
  44     {
  45       if(!fParamSpace[i]) continue;
  46       delete fParamSpace[i];
  47     }
  48   delete [] fParamSpace;
  49 }
  50
  51 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
  52                                              Int_t nybin,Double_t phimin,Double_t phimax)
  53 {
  54   //Create the histograms (parameter space).
  55   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  56   //The arguments give the range and binning;
  57   //nxbin = #bins in kappa
  58   //nybin = #bins in phi0
  59   //pt_min = mimium Pt of track (corresponding to maximum kappa)
  60   //phi_min = mimimum phi0 (degrees)
  61   //phi_max = maximum phi0 (degrees)
  62
  63   Double_t bfact = 0.0029980;
  64   Double_t bfield = 0.2;
  65   Double_t x = bfact*bfield/pt_min;
  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   if(!fTransform)
 130     {
 131       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Transformer")
 132         <<"No AliL3Transform object"<<ENDLOG;
 133       return;
 134     }
 135
 136   //Loop over the padrows:
 137   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
 138     {
 139       //Get the data on this padrow:
 140       AliL3DigitData *digPt = tempPt->fDigitData;
 141       if(i != (Int_t)tempPt->fRow)
 142         {
 143           printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
 144           continue;
 145         }
 146
 147       //Loop over the data on this padrow:
 148       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 149         {
 150           UShort_t charge = digPt[j].fCharge;
 151           UChar_t pad = digPt[j].fPad;
 152           UShort_t time = digPt[j].fTime;
 153           if(charge <= GetThreshold())
 154             continue;
 155           Int_t sector,row;
 156           Float_t xyz[3];
 157
 158           //Transform data to local cartesian coordinates:
 159           fTransform->Slice2Sector(GetSlice(),i,sector,row);
 160           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 161
 162           //Calculate the eta:
 163           Double_t eta = fTransform->GetEta(xyz);
 164
 165           //Get the corresponding index, which determines which histogram to fill:
 166           Int_t eta_index = GetEtaIndex(eta);
 167           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 168             continue;
 169
 170           //Get the correct histogrampointer:
 171           AliL3Histogram *hist = fParamSpace[eta_index];
 172           if(!hist)
 173             {
 174               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
 175               continue;
 176             }
 177
 178           //Do the transformation:
 179           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 180           Float_t phi = fTransform->GetPhi(xyz);
 181
 182           //Fill the histogram along the phirange
 183           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 184             {
 185               Float_t phi0 = hist->GetBinCenterY(b);
 186               Float_t kappa = 2*sin(phi - phi0)/R;
 187               hist->Fill(kappa,phi0,charge);
 188             }
 189         }
 190
 191       //Move the data pointer to the next padrow:
 192       AliL3MemHandler::UpdateRowPointer(tempPt);
 193     }
 194 }
 195
 196 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
 197 {
 198   //Circle transform, using combinations of every 2 points lying
 199   //on different padrows and within the same etaslice.
 200
 201   AliL3DigitRowData *tempPt = GetDataPointer();
 202   if(!tempPt)
 203     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 204       <<"No input data "<<ENDLOG;
 205
 206   if(!fTransform)
 207     {
 208       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Transformer")
 209         <<"No AliL3Transform object"<<ENDLOG;
 210       return;
 211     }
 212
 213   Int_t counter=0;
 214   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
 215     {
 216       counter += tempPt->fNDigit;
 217       AliL3MemHandler::UpdateRowPointer(tempPt);
 218     }
 219
 220   struct Digit {
 221     Int_t row;
 222     Double_t r;
 223     Double_t phi;
 224     Int_t eta_index;
 225     Int_t charge;
 226   };
 227
 228   Digit *digits = new Digit[counter];
 229   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
 230
 231   Int_t total_digits=counter;
 232   Int_t sector,row,tot_charge,pad,time,charge;
 233   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
 234   Float_t xyz[3];
 235
 236   counter=0;
 237   tempPt = GetDataPointer();
 238
 239   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
 240     {
 241       AliL3DigitData *digPt = tempPt->fDigitData;
 242       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 243         {
 244           charge = digPt[di].fCharge;
 245           pad = digPt[di].fPad;
 246           time = digPt[di].fTime;
 247           fTransform->Slice2Sector(GetSlice(),i,sector,row);
 248           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 249           eta = fTransform->GetEta(xyz);
 250           digits[counter].row = i;
 251           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 252           digits[counter].phi = atan2(xyz[1],xyz[0]);
 253           digits[counter].eta_index = GetEtaIndex(eta);
 254           digits[counter].charge = charge;
 255           counter++;
 256         }
 257       AliL3MemHandler::UpdateRowPointer(tempPt);
 258     }
 259
 260   for(Int_t i=0; i<total_digits; i++)
 261     {
 262       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
 263       Int_t ind = digits[i].eta_index;
 264
 265       for(Int_t j=i+1; j<total_digits; j++)
 266         {
 267           if(digits[i].row == digits[j].row) continue;
 268           if(digits[i].eta_index != digits[j].eta_index) continue;
 269           if(digits[i].row + row_range < digits[j].row) break;
 270
 271           //Get the correct histogrampointer:
 272           AliL3Histogram *hist = fParamSpace[ind];
 273           if(!hist)
 274             {
 275               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
 276               continue;
 277             }
 278
 279           r1 = digits[i].r;
 280           phi1 = digits[i].phi;
 281           r2 = digits[j].r;
 282           phi2 = digits[j].phi;
 283           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
 284           kappa = 2*sin(phi2-phi_0)/r2;
 285           tot_charge = digits[i].charge + digits[j].charge;
 286           hist->Fill(kappa,phi_0,tot_charge);
 287         }
 288     }
 289   delete [] digits;
 290 }
 291
 292 void AliL3HoughTransformer::TransformLine()
 293 {
 294   //Do a line transform on the data.
 295
 296
 297   AliL3DigitRowData *tempPt = GetDataPointer();
 298   if(!tempPt)
 299     {
 300       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
 301         <<"No input data "<<ENDLOG;
 302       return;
 303     }
 304   if(!fTransform)
 305     {
 306       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Transformer")
 307         <<"No AliL3Transform object"<<ENDLOG;
 308       return;
 309     }
 310
 311   for(Int_t i=NRows[GetPatch()][0]; i<=NRows[GetPatch()][1]; i++)
 312     {
 313       AliL3DigitData *digPt = tempPt->fDigitData;
 314       if(i != (Int_t)tempPt->fRow)
 315         {
 316           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
 317           continue;
 318         }
 319       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 320         {
 321           UShort_t charge = digPt[j].fCharge;
 322           UChar_t pad = digPt[j].fPad;
 323           UShort_t time = digPt[j].fTime;
 324           if(charge < GetThreshold())
 325             continue;
 326           Int_t sector,row;
 327           Float_t xyz[3];
 328           fTransform->Slice2Sector(GetSlice(),i,sector,row);
 329           fTransform->Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 330           Float_t eta = fTransform->GetEta(xyz);
 331           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 332           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 333             continue;
 334
 335           //Get the correct histogram:
 336           AliL3Histogram *hist = fParamSpace[eta_index];
 337           if(!hist)
 338             {
 339               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
 340               continue;
 341             }
 342           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
 343             {
 344               Double_t theta = hist->GetBinCenterX(xbin);
 345               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 346               hist->Fill(theta,rho,charge);
 347             }
 348         }
 349       AliL3MemHandler::UpdateRowPointer(tempPt);
 350     }
 351
 352 }