HLT/hough/AliL3HoughClusterTransformer.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Logging.h"
   9 #include "AliL3HoughClusterTransformer.h"
  10 #include "AliL3MemHandler.h"
  11 #include "AliL3SpacePointData.h"
  12 #include "AliL3Transform.h"
  13 #include "AliL3DigitData.h"
  14 #include "AliL3Histogram.h"
  15 #include "AliL3ClustFinderNew.h"
  16
  17 #if __GNUC__ == 3
  18 using namespace std;
  19 #endif
  20
  21 //_____________________________________________________________
  22 // AliL3HoughClusterTransformer
  23 //
  24 // Hough transformation class.
  25 //
  26 // This class performs the Hough transform on _clusters_, and not raw data.
  27 // Therefore, it first finds the clusters using the HLT cluster finder, and
  28 // then uses these as input for the transform.
  29
  30 ClassImp(AliL3HoughClusterTransformer)
  31
  32 AliL3HoughClusterTransformer::AliL3HoughClusterTransformer()
  33 {
  34   //Default constructor
  35   fParamSpace = 0;
  36   fNClusters = 0;
  37   fClusters = 0;
  38   fMemHandler = 0;
  39 #ifdef do_mc
  40   fTrackID = 0;
  41 #endif
  42 }
  43
  44 AliL3HoughClusterTransformer::AliL3HoughClusterTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
  45 {
  46   //Normal constructor
  47   fParamSpace = 0;
  48   fNClusters=0;
  49   fClusters=0;
  50   fMemHandler=0;
  51 #ifdef do_mc
  52   fTrackID=0;
  53 #endif
  54 }
  55
  56 AliL3HoughClusterTransformer::~AliL3HoughClusterTransformer()
  57 {
  58   DeleteHistograms();
  59   if(fMemHandler)
  60     delete fMemHandler;
  61 #ifdef do_mc
  62   if(fTrackID)
  63     {
  64       for(Int_t i=0; i<GetNEtaSegments(); i++)
  65         {
  66           if(!fTrackID[i]) continue;
  67           delete fTrackID[i];
  68         }
  69       delete [] fTrackID;
  70     }
  71 #endif
  72 }
  73
  74 //void AliL3HoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
  75
  76 void AliL3HoughClusterTransformer::DeleteHistograms()
  77 {
  78   if(fParamSpace)
  79     {
  80       for(Int_t i=0; i<GetNEtaSegments(); i++)
  81         {
  82           if(!fParamSpace[i]) continue;
  83           delete fParamSpace[i];
  84         }
  85       delete [] fParamSpace;
  86     }
  87 }
  88
  89 void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t pt_min,
  90                                              Int_t nybin,Float_t phimin,Float_t phimax)
  91 {
  92   //Create the histograms (parameter space).
  93   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  94   //The arguments give the range and binning;
  95   //nxbin = #bins in kappa
  96   //nybin = #bins in phi0
  97   //pt_min = mimium Pt of track (corresponding to maximum kappa)
  98   //phi_min = mimimum phi0 (degrees)
  99   //phi_max = maximum phi0 (degrees)
 100
 101   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
 102   Double_t torad = AliL3Transform::Pi()/180;
 103   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
 104 }
 105
 106 void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
 107                                                     Int_t nybin,Float_t ymin,Float_t ymax)
 108 {
 109   Char_t histname[256];
 110
 111   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
 112   for(Int_t i=0; i<GetNEtaSegments(); i++)
 113     {
 114       sprintf(histname,"paramspace_%d",i);
 115       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 116     }
 117 #ifdef do_mc
 118   Int_t ncells = (nxbin+2)*(nybin+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 #endif
 124
 125 }
 126
 127 void AliL3HoughClusterTransformer::Reset()
 128 {
 129   //Reset all the histograms
 130
 131   if(!fParamSpace)
 132     {
 133       LOG(AliL3Log::kWarning,"AliL3HoughClusterTransformer::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
 141   if(fMemHandler)
 142     delete fMemHandler;
 143   fNClusters=0;
 144   fClusters=0;
 145 #ifdef do_mc
 146   AliL3Histogram *hist = fParamSpace[0];
 147   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 148   for(Int_t i=0; i<GetNEtaSegments(); i++)
 149     memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
 150 #endif
 151 }
 152
 153 Int_t AliL3HoughClusterTransformer::GetEtaIndex(Double_t eta)
 154 {
 155   //Return the histogram index of the corresponding eta.
 156
 157   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 158   Double_t index = (eta-GetEtaMin())/etaslice;
 159   return (Int_t)index;
 160 }
 161
 162 inline AliL3Histogram *AliL3HoughClusterTransformer::GetHistogram(Int_t eta_index)
 163 {
 164   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
 165     return 0;
 166   if(!fParamSpace[eta_index])
 167     return 0;
 168   return fParamSpace[eta_index];
 169 }
 170
 171 Double_t AliL3HoughClusterTransformer::GetEta(Int_t eta_index,Int_t slice)
 172 {
 173   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 174   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
 175   if(slice>17) eta*=-1;
 176   return eta;
 177 }
 178
 179 void AliL3HoughClusterTransformer::FindClusters()
 180 {
 181   //Find the clusters
 182   if(!GetDataPointer())
 183     {
 184       cerr<<"AliL3HoughClusterTransformer::FindClusters : Zero data pointer"<<endl;
 185       return;
 186     }
 187
 188   const Int_t maxpoints=100000;
 189   const Int_t pointsize = maxpoints * sizeof(AliL3SpacePointData);
 190
 191   fMemHandler = new AliL3MemHandler();
 192   fClusters = (AliL3SpacePointData*)fMemHandler->Allocate(pointsize);
 193   AliL3ClustFinderNew *cf = new AliL3ClustFinderNew();
 194   cf->InitSlice(0,GetPatch(),AliL3Transform::GetFirstRow(GetPatch()),AliL3Transform::GetLastRow(GetPatch()),maxpoints);
 195   cf->SetDeconv(kFALSE);
 196   cf->SetOutputArray(fClusters);
 197   cf->Read(1,GetDataPointer());
 198   cf->ProcessDigits();
 199   fNClusters = cf->GetNumberOfClusters();
 200   delete cf;
 201 }
 202
 203 void AliL3HoughClusterTransformer::TransformCircle()
 204 {
 205   //Transform the input data with a circle HT.
 206   //The function loops over all the data, and transforms each cluster with the equations:
 207   //
 208   //kappa = 2/R*sin(phi - phi0)
 209   //
 210   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
 211   //
 212   //Each cluster then transforms into a curve in the (kappa,phi0)-space. In order to find
 213   //which histogram in which the pixel should be transformed, the eta-value is calcluated
 214   //and the proper histogram index is found by GetEtaIndex(eta).
 215
 216   FindClusters();
 217   if(!fClusters)
 218     {
 219       LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircle","Data")
 220         <<"No input data "<<ENDLOG;
 221       return;
 222     }
 223
 224   Float_t xyz[3];
 225   for(Int_t i=0; i<fNClusters; i++)
 226     {
 227       xyz[0] = fClusters[i].fX;
 228       xyz[1] = fClusters[i].fY;
 229       xyz[2] = fClusters[i].fZ;
 230       Double_t eta = AliL3Transform::GetEta(xyz);
 231       Int_t eta_index = GetEtaIndex(eta);
 232       if(eta_index < 0 || eta_index >= GetNEtaSegments())
 233         continue;
 234
 235       AliL3Histogram *hist = fParamSpace[eta_index];
 236
 237       //Do the transformation:
 238       Double_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 239       Double_t phi = AliL3Transform::GetPhi(xyz);
 240
 241       //Fill the histogram along the phirange
 242       for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 243         {
 244           Double_t phi0 = hist->GetBinCenterY(b);
 245           Double_t kappa = 2*sin(phi - phi0)/R;
 246           hist->Fill(kappa,phi0,1);
 247 #ifdef do_mc
 248           Int_t bin = hist->FindBin(kappa,phi0);
 249           for(Int_t t=0; t<3; t++)
 250             {
 251               Int_t label = fClusters[i].fTrackID[t];
 252               if(label < 0) break;
 253               UInt_t c;
 254               for(c=0; c<MaxTrack; c++)
 255                 if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
 256                   break;
 257               if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 258               fTrackID[eta_index][bin].fLabel[c] = label;
 259               fTrackID[eta_index][bin].fNHits[c]++;
 260             }
 261 #endif
 262         }
 263     }
 264 }
 265
 266 void AliL3HoughClusterTransformer::TransformCircleC(Int_t *row_range,Int_t every)
 267 {
 268   //Circle transform, using combinations of every 2 points lying
 269   //on different padrows and within the same etaslice.
 270
 271   FindClusters();
 272   if(!fClusters)
 273     LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircleC","Data")
 274       <<"No input data "<<ENDLOG;
 275
 276   Float_t xyz[3];
 277   Double_t eta,r1,phi1,r2,phi2,kappa,psi;
 278   Int_t index1,index2;
 279   for(Int_t i=0; i<fNClusters; i++)
 280     {
 281       xyz[0] = fClusters[i].fX;
 282       xyz[1] = fClusters[i].fY;
 283       xyz[2] = fClusters[i].fZ;
 284       r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 285       phi1 = atan2(xyz[1],xyz[0]);
 286       eta = AliL3Transform::GetEta(xyz);
 287       index1 = GetEtaIndex(eta);
 288       if(index1 < 0 || index1 >= GetNEtaSegments()) continue;
 289
 290       for(Int_t j=i+1; j<fNClusters; j++)
 291         {
 292           if(fClusters[j].fPadRow == fClusters[i].fPadRow) continue;
 293           xyz[0] = fClusters[j].fX;
 294           xyz[1] = fClusters[j].fY;
 295           xyz[2] = fClusters[j].fZ;
 296           eta = AliL3Transform::GetEta(xyz);
 297           index2 = GetEtaIndex(eta);
 298           if(index1 != index2) continue;
 299
 300           AliL3Histogram *hist = fParamSpace[index1];
 301
 302           r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 303           phi2 = atan2(xyz[1],xyz[0]);
 304
 305           psi = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
 306           kappa = 2*sin(phi1-psi) / r1;
 307           hist->Fill(kappa,psi,1);
 308 #ifdef do_mc
 309           Int_t bin = hist->FindBin(kappa,psi);
 310           for(Int_t l=0; l<3; l++)
 311             {
 312               for(Int_t m=0; m<3; m++)
 313                 {
 314                   if(fClusters[i].fTrackID[l] == fClusters[j].fTrackID[m])
 315                     {
 316                       Int_t label = fClusters[i].fTrackID[l];
 317                       if(label < 0) continue;
 318                       UInt_t c;
 319                       for(c=0; c<MaxTrack; c++)
 320                         if(fTrackID[index1][bin].fLabel[c] == label || fTrackID[index1][bin].fNHits[c] == 0)
 321                           break;
 322                       if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
 323                       fTrackID[index1][bin].fLabel[c] = label;
 324                       fTrackID[index1][bin].fNHits[c]++;
 325                     }
 326                 }
 327             }
 328 #endif
 329         }
 330     }
 331 }
 332
 333
 334 Int_t AliL3HoughClusterTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
 335 {
 336
 337 #ifdef do_mc
 338   if(eta_index < 0 || eta_index > GetNEtaSegments())
 339     {
 340       cerr<<"AliL3HoughClusterTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
 341       return -1;
 342     }
 343   AliL3Histogram *hist = fParamSpace[eta_index];
 344   Int_t bin = hist->FindBin(kappa,psi);
 345   Int_t label=-1;
 346   Int_t max=0;
 347   for(UInt_t i=0; i<MaxTrack; i++)
 348     {
 349       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
 350       if(nhits == 0) break;
 351       if(nhits > max)
 352         {
 353           max = nhits;
 354           label = fTrackID[eta_index][bin].fLabel[i];
 355         }
 356     }
 357   return label;
 358 #endif
 359   cout<<"AliL3HoughClusterTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
 360   return -1;
 361 }
 362