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 netasegments) : AliL3HoughBaseTransformer(slice,patch,netasegments)
  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   //dtor
  59   DeleteHistograms();
  60   if(fMemHandler)
  61     delete fMemHandler;
  62 #ifdef do_mc
  63   if(fTrackID)
  64     {
  65       for(Int_t i=0; i<GetNEtaSegments(); i++)
  66         {
  67           if(!fTrackID[i]) continue;
  68           delete fTrackID[i];
  69         }
  70       delete [] fTrackID;
  71     }
  72 #endif
  73 }
  74
  75 //void AliL3HoughClusterTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t netasegments=100){}
  76
  77 void AliL3HoughClusterTransformer::DeleteHistograms()
  78 {
  79   //Delete hough space histograms
  80   if(fParamSpace)
  81     {
  82       for(Int_t i=0; i<GetNEtaSegments(); i++)
  83         {
  84           if(!fParamSpace[i]) continue;
  85           delete fParamSpace[i];
  86         }
  87       delete [] fParamSpace;
  88     }
  89 }
  90
  91 void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t ptmin,
  92                                              Int_t nybin,Float_t phimin,Float_t phimax)
  93 {
  94   //Create the histograms (parameter space).
  95   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  96   //The arguments give the range and binning;
  97   //nxbin = #bins in kappa
  98   //nybin = #bins in phi0
  99   //ptmin = mimium Pt of track (corresponding to maximum kappa)
 100   //phimin = mimimum phi0 (degrees)
 101   //phimax = maximum phi0 (degrees)
 102
 103   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/ptmin;
 104   Double_t torad = AliL3Transform::Pi()/180;
 105   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
 106 }
 107
 108 void AliL3HoughClusterTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
 109                                                     Int_t nybin,Float_t ymin,Float_t ymax)
 110 {
 111   //Create histograms which contain hough space
 112   Char_t histname[256];
 113
 114   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
 115   for(Int_t i=0; i<GetNEtaSegments(); i++)
 116     {
 117       sprintf(histname,"paramspace_%d",i);
 118       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 119     }
 120 #ifdef do_mc
 121   Int_t ncells = (nxbin+2)*(nybin+2);
 122   cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fTrackID"<<endl;
 123   fTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
 124   for(Int_t i=0; i<GetNEtaSegments(); i++)
 125     fTrackID[i] = new AliL3TrackIndex[ncells];
 126 #endif
 127
 128 }
 129
 130 void AliL3HoughClusterTransformer::Reset()
 131 {
 132   //Reset all the histograms
 133
 134   if(!fParamSpace)
 135     {
 136       LOG(AliL3Log::kWarning,"AliL3HoughClusterTransformer::Reset","Histograms")
 137         <<"No histograms to reset"<<ENDLOG;
 138       return;
 139     }
 140
 141   for(Int_t i=0; i<GetNEtaSegments(); i++)
 142     fParamSpace[i]->Reset();
 143
 144   if(fMemHandler)
 145     delete fMemHandler;
 146   fNClusters=0;
 147   fClusters=0;
 148 #ifdef do_mc
 149   AliL3Histogram *hist = fParamSpace[0];
 150   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 151   for(Int_t i=0; i<GetNEtaSegments(); i++)
 152     memset(fTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
 153 #endif
 154 }
 155
 156 Int_t AliL3HoughClusterTransformer::GetEtaIndex(Double_t eta) const
 157 {
 158   //Return the histogram index of the corresponding eta.
 159
 160   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 161   Double_t index = (eta-GetEtaMin())/etaslice;
 162   return (Int_t)index;
 163 }
 164
 165 inline AliL3Histogram *AliL3HoughClusterTransformer::GetHistogram(Int_t etaindex)
 166 {
 167   //Returns the histogram which correspond to a given eta slice
 168   if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
 169     return 0;
 170   if(!fParamSpace[etaindex])
 171     return 0;
 172   return fParamSpace[etaindex];
 173 }
 174
 175 Double_t AliL3HoughClusterTransformer::GetEta(Int_t etaindex,Int_t slice) const
 176 {
 177   //Returns eta associated with given eta slice
 178   Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 179   Double_t eta=(Double_t)((etaindex+0.5)*etaslice);
 180   if(slice>17) eta*=-1;
 181   return eta;
 182 }
 183
 184 void AliL3HoughClusterTransformer::FindClusters()
 185 {
 186   //Find the clusters
 187   if(!GetDataPointer())
 188     {
 189       cerr<<"AliL3HoughClusterTransformer::FindClusters : Zero data pointer"<<endl;
 190       return;
 191     }
 192
 193   const Int_t kMaxpoints=100000;
 194   const Int_t kPointsize = kMaxpoints * sizeof(AliL3SpacePointData);
 195
 196   fMemHandler = new AliL3MemHandler();
 197   fClusters = (AliL3SpacePointData*)fMemHandler->Allocate(kPointsize);
 198   AliL3ClustFinderNew *cf = new AliL3ClustFinderNew();
 199   cf->InitSlice(0,GetPatch(),AliL3Transform::GetFirstRow(GetPatch()),AliL3Transform::GetLastRow(GetPatch()),kMaxpoints);
 200   cf->SetDeconv(kFALSE);
 201   cf->SetOutputArray(fClusters);
 202   cf->Read(1,GetDataPointer());
 203   cf->ProcessDigits();
 204   fNClusters = cf->GetNumberOfClusters();
 205   delete cf;
 206 }
 207
 208 void AliL3HoughClusterTransformer::TransformCircle()
 209 {
 210   //Transform the input data with a circle HT.
 211   //The function loops over all the data, and transforms each cluster with the equations:
 212   //
 213   //kappa = 2/r*sin(phi - phi0)
 214   //
 215   //where r = sqrt(x*x +y*y), and phi = arctan(y/x)
 216   //
 217   //Each cluster then transforms into a curve in the (kappa,phi0)-space. In order to find
 218   //which histogram in which the pixel should be transformed, the eta-value is calcluated
 219   //and the proper histogram index is found by GetEtaIndex(eta).
 220
 221   FindClusters();
 222   if(!fClusters)
 223     {
 224       LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircle","Data")
 225         <<"No input data "<<ENDLOG;
 226       return;
 227     }
 228
 229   Float_t xyz[3];
 230   for(Int_t i=0; i<fNClusters; i++)
 231     {
 232       xyz[0] = fClusters[i].fX;
 233       xyz[1] = fClusters[i].fY;
 234       xyz[2] = fClusters[i].fZ;
 235       Double_t eta = AliL3Transform::GetEta(xyz);
 236       Int_t etaindex = GetEtaIndex(eta);
 237       if(etaindex < 0 || etaindex >= GetNEtaSegments())
 238         continue;
 239
 240       AliL3Histogram *hist = fParamSpace[etaindex];
 241
 242       //Do the transformation:
 243       Double_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 244       Double_t phi = AliL3Transform::GetPhi(xyz);
 245
 246       //Fill the histogram along the phirange
 247       for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 248         {
 249           Double_t phi0 = hist->GetBinCenterY(b);
 250           Double_t kappa = 2*sin(phi - phi0)/r;
 251           hist->Fill(kappa,phi0,1);
 252 #ifdef do_mc
 253           Int_t bin = hist->FindBin(kappa,phi0);
 254           for(Int_t t=0; t<3; t++)
 255             {
 256               Int_t label = fClusters[i].fTrackID[t];
 257               if(label < 0) break;
 258               UInt_t c;
 259               for(c=0; c<MaxTrack; c++)
 260                 if(fTrackID[etaindex][bin].fLabel[c] == label || fTrackID[etaindex][bin].fNHits[c] == 0)
 261                   break;
 262               if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 263               fTrackID[etaindex][bin].fLabel[c] = label;
 264               fTrackID[etaindex][bin].fNHits[c]++;
 265             }
 266 #endif
 267         }
 268     }
 269 }
 270
 271 void AliL3HoughClusterTransformer::TransformCircleC(Int_t */*row_range*/,Int_t /*every*/)
 272 {
 273   //Circle transform, using combinations of every 2 points lying
 274   //on different padrows and within the same etaslice.
 275
 276   FindClusters();
 277   if(!fClusters)
 278     LOG(AliL3Log::kError,"AliL3HoughClusterTransformer::TransformCircleC","Data")
 279       <<"No input data "<<ENDLOG;
 280
 281   Float_t xyz[3];
 282   Double_t eta,r1,phi1,r2,phi2,kappa,psi;
 283   Int_t index1,index2;
 284   for(Int_t i=0; i<fNClusters; i++)
 285     {
 286       xyz[0] = fClusters[i].fX;
 287       xyz[1] = fClusters[i].fY;
 288       xyz[2] = fClusters[i].fZ;
 289       r1 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 290       phi1 = atan2(xyz[1],xyz[0]);
 291       eta = AliL3Transform::GetEta(xyz);
 292       index1 = GetEtaIndex(eta);
 293       if(index1 < 0 || index1 >= GetNEtaSegments()) continue;
 294
 295       for(Int_t j=i+1; j<fNClusters; j++)
 296         {
 297           if(fClusters[j].fPadRow == fClusters[i].fPadRow) continue;
 298           xyz[0] = fClusters[j].fX;
 299           xyz[1] = fClusters[j].fY;
 300           xyz[2] = fClusters[j].fZ;
 301           eta = AliL3Transform::GetEta(xyz);
 302           index2 = GetEtaIndex(eta);
 303           if(index1 != index2) continue;
 304
 305           AliL3Histogram *hist = fParamSpace[index1];
 306
 307           r2 = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 308           phi2 = atan2(xyz[1],xyz[0]);
 309
 310           psi = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
 311           kappa = 2*sin(phi1-psi) / r1;
 312           hist->Fill(kappa,psi,1);
 313 #ifdef do_mc
 314           Int_t bin = hist->FindBin(kappa,psi);
 315           for(Int_t l=0; l<3; l++)
 316             {
 317               for(Int_t m=0; m<3; m++)
 318                 {
 319                   if(fClusters[i].fTrackID[l] == fClusters[j].fTrackID[m])
 320                     {
 321                       Int_t label = fClusters[i].fTrackID[l];
 322                       if(label < 0) continue;
 323                       UInt_t c;
 324                       for(c=0; c<MaxTrack; c++)
 325                         if(fTrackID[index1][bin].fLabel[c] == label || fTrackID[index1][bin].fNHits[c] == 0)
 326                           break;
 327                       if(c == MaxTrack-1) cerr<<"AliL3HoughClusterTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
 328                       fTrackID[index1][bin].fLabel[c] = label;
 329                       fTrackID[index1][bin].fNHits[c]++;
 330                     }
 331                 }
 332             }
 333 #endif
 334         }
 335     }
 336 }
 337
 338
 339 #ifdef do_mc
 340 Int_t AliL3HoughClusterTransformer::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi) const
 341 {
 342   //Returns the mc label for a given bin in the hough space
 343   if(etaindex < 0 || etaindex > GetNEtaSegments())
 344     {
 345       cerr<<"AliL3HoughClusterTransformer::GetTrackID : Wrong etaindex "<<etaindex<<endl;
 346       return -1;
 347     }
 348   AliL3Histogram *hist = fParamSpace[etaindex];
 349   Int_t bin = hist->FindBin(kappa,psi);
 350   Int_t label=-1;
 351   Int_t max=0;
 352   for(UInt_t i=0; i<MaxTrack; i++)
 353     {
 354       Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
 355       if(nhits == 0) break;
 356       if(nhits > max)
 357         {
 358           max = nhits;
 359           label = fTrackID[etaindex][bin].fLabel[i];
 360         }
 361     }
 362   return label;
 363 #else
 364   Int_t AliL3HoughClusterTransformer::GetTrackID(Int_t /*etaindex*/,Double_t /*kappa*/,Double_t /*psi*/) const
 365 {
 366   // Does nothing if do_mc undefinde
 367   cout<<"AliL3HoughClusterTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
 368   return -1;
 369 #endif
 370 }
 371