TRD/AliTRDchamberTimeBin.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id: AliTRDchamberTimeBin.cxx 23313 2008-01-11 14:56:43Z cblume $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Organization of clusters at the level of 1 TRD chamber.                  //
  21 //  The data structure is used for tracking at the stack level.              //
  22 //                                                                           //
  23 //  Functionalities:                                                         //
  24 //   1. cluster organization and sorting                                     //
  25 //   2. fast data navigation                                                 //
  26 //                                                                           //
  27 //  Authors:                                                                 //
  28 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  29 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  30 //                                                                           //
  31 ///////////////////////////////////////////////////////////////////////////////
  32
  33 #include <TObject.h>
  34 #include <TROOT.h>
  35 #include <TMath.h>
  36 #include <TStopwatch.h>
  37 #include <TTreeStream.h>
  38
  39 #include "AliLog.h"
  40
  41 #include "AliTRDgeometry.h"
  42 #include "AliTRDpadPlane.h"
  43 #include "AliTRDchamberTimeBin.h"
  44 #include "AliTRDrecoParam.h"
  45 #include "AliTRDReconstructor.h"
  46
  47 ClassImp(AliTRDchamberTimeBin)
  48
  49 //_____________________________________________________________________________
  50 AliTRDchamberTimeBin::AliTRDchamberTimeBin(Int_t plane, Int_t stack, Int_t sector, Double_t z0, Double_t zLength)
  51   :TObject()
  52   ,fReconstructor(0x0)
  53   ,fPlane(plane)
  54   ,fStack(stack)
  55   ,fSector(sector)
  56   ,fNRows(kMaxRows)
  57   ,fN(0)
  58   ,fX(0.)
  59   ,fZ0(z0)
  60   ,fZLength(zLength)
  61 {
  62   //
  63   // Default constructor (Only provided to use AliTRDchamberTimeBin with arrays)
  64   //
  65   SetBit(kT0, kFALSE);
  66   SetBit(kOwner, kFALSE);
  67   memset(fPositions, 1, kMaxRows*sizeof(UChar_t));
  68   memset(fClusters, 0, kMaxClustersLayer*sizeof(AliTRDcluster*));
  69   memset(fIndex, 1, kMaxClustersLayer*sizeof(UInt_t));
  70 }
  71
  72
  73 //_____________________________________________________________________________
  74 AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer):
  75   TObject()
  76   ,fReconstructor(layer.fReconstructor)
  77   ,fPlane(layer.fPlane)
  78   ,fStack(layer.fStack)
  79   ,fSector(layer.fSector)
  80   ,fNRows(layer.fNRows)
  81   ,fN(layer.fN)
  82   ,fX(layer.fX)
  83   ,fZ0(layer.fZ0)
  84   ,fZLength(layer.fZLength)
  85 {
  86 // Copy Constructor
  87
  88   SetBit(kT0, layer.IsT0());
  89   SetBit(kOwner, kFALSE);
  90   for(int i=0; i<kMaxRows; i++) fPositions[i] = layer.fPositions[i];
  91   memcpy(&fClusters[0], &layer.fClusters[0], kMaxClustersLayer*sizeof(AliTRDcluster*));
  92   memcpy(&fIndex[0], &layer.fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
  93
  94
  95 //      BuildIndices();
  96 }
  97
  98 //_____________________________________________________________________________
  99 AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDchamberTimeBin &layer)
 100 {
 101 // Assignment operator
 102
 103   if (this != &layer) layer.Copy(*this);
 104   return *this;
 105 }
 106
 107 //_____________________________________________________________________________
 108 void AliTRDchamberTimeBin::Clear(const Option_t *)
 109 {
 110   for(Int_t it = 0; it<kMaxClustersLayer; it++){
 111     if(IsOwner()) delete fClusters[it];
 112     fClusters[it] = NULL;
 113   }
 114   fN = 0;
 115 }
 116
 117 //_____________________________________________________________________________
 118 void AliTRDchamberTimeBin::Copy(TObject &o) const
 119 {
 120 // Copy method. Performs a deep copy of all data from this object to object o.
 121
 122   AliTRDchamberTimeBin &layer = (AliTRDchamberTimeBin &)o;
 123   layer.fReconstructor = fReconstructor;
 124   layer.fPlane       = fPlane;
 125   layer.fStack       = fStack;
 126   layer.fSector      = fSector;
 127   layer.fNRows       = fNRows;
 128   layer.fN           = fN;
 129   layer.fX           = fX;
 130   layer.fZ0          = fZ0;
 131   layer.fZLength     = fZLength;
 132   layer.SetT0(IsT0());
 133   layer.SetBit(kOwner, kFALSE);
 134
 135   for(int i=0; i<kMaxRows; i++) layer.fPositions[i] = fPositions[i];
 136   memcpy(&layer.fClusters[0], &fClusters[0], kMaxClustersLayer*sizeof(AliTRDcluster*));
 137   memcpy(&layer.fIndex[0], &fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
 138
 139   TObject::Copy(layer); // copies everything into layer
 140
 141 //      layer.BuildIndices();
 142 }
 143
 144 //_____________________________________________________________________________
 145 AliTRDchamberTimeBin::~AliTRDchamberTimeBin()
 146 {
 147 // Destructor
 148   if(IsOwner()){
 149     for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++) delete (*cit);
 150   }
 151 }
 152
 153 //_____________________________________________________________________________
 154 void  AliTRDchamberTimeBin::SetOwner(Bool_t copy)
 155 {
 156 // Sets the ownership of the clusters to this
 157 // If option "copy" is kTRUE [default] the clusters
 158 // are also copied otherwise only the ownership bit
 159 // is flipped.
 160
 161   SetBit(kOwner, kTRUE);
 162   if(!copy) return;
 163   for(AliTRDcluster **cit = &fClusters[0]; (*cit); ++cit){
 164     (*cit) = new AliTRDcluster(*(*cit));
 165   }
 166 }
 167
 168 //_____________________________________________________________________________
 169 void AliTRDchamberTimeBin::SetRange(Float_t z0, Float_t zLength)
 170 {
 171 // Sets the range in z-direction
 172 //
 173 // Parameters:
 174 //   z0      : starting position of layer in the z direction
 175 //   zLength : length of layer in the z direction
 176
 177   fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
 178   fZLength = TMath::Abs(zLength);
 179 }
 180
 181 //_____________________________________________________________________________
 182 void AliTRDchamberTimeBin::InsertCluster(AliTRDcluster *c, UInt_t index)
 183 {
 184   //
 185   // Insert cluster in cluster array.
 186   // Clusters are sorted according to Y coordinate.
 187   //
 188
 189   //if (fTimeBinIndex < 0) {
 190     //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
 191     //return;
 192   //}
 193
 194   if (fN == (Int_t) kMaxClustersLayer) {
 195     //AliWarning("Too many clusters !\n");
 196     return;
 197   }
 198
 199   if (fN == 0) {
 200     fIndex[0]       = index;
 201     fClusters[fN++] = c;
 202     return;
 203   }
 204
 205   Int_t i = Find(c->GetY());
 206   memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
 207   memmove(fIndex   +i+1,fIndex   +i,(fN-i)*sizeof(UInt_t));
 208   fIndex[i]    = index;
 209   fClusters[i] = c;
 210   fN++;
 211 }
 212
 213 //___________________________________________________
 214 void AliTRDchamberTimeBin::Bootstrap(const AliTRDReconstructor *rec, Int_t det)
 215 {
 216 // Reinitialize all data members from the clusters array
 217 // It has to be used after reading from disk
 218
 219   fReconstructor = rec;
 220   fPlane  = AliTRDgeometry::GetLayer(det);
 221   fStack  = AliTRDgeometry::GetStack(det);
 222   fSector = AliTRDgeometry::GetSector(det);
 223   AliTRDgeometry g;
 224   fNRows  = g.GetPadPlane(fPlane, fStack)->GetNrows();
 225   fN = 0;
 226   for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++) fN++;
 227   BuildIndices();
 228 }
 229
 230 //_____________________________________________________________________________
 231 void AliTRDchamberTimeBin::BuildIndices(Int_t iter)
 232 {
 233 // Rearrangement of the clusters belonging to the propagation layer for the stack.
 234 //
 235 // Detailed description
 236 //
 237 // The array indices of all clusters in one PropagationLayer are stored in
 238 // array. The array is divided into several bins.
 239 // The clusters are sorted in increasing order of their y coordinate.
 240 //
 241 // Sorting algorithm: TreeSearch
 242 //
 243
 244   if(!fN) return;
 245
 246   // Select clusters that belong to the Stack
 247   Int_t nClStack = 0;                                   // Internal counter
 248   for(Int_t i = 0; i < fN; i++){
 249     if(fClusters[i]->IsUsed() || fClusters[i]->IsShared()){
 250       fClusters[i] = 0x0;
 251       fIndex[i] = 0xffff;
 252     } else nClStack++;
 253   }
 254   if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d. Truncating.", nClStack, kMaxClustersLayer));
 255
 256   // Nothing in this time bin. Reset indexes
 257   if(!nClStack){
 258     fN = 0;
 259     memset(&fPositions[0], 0xff, sizeof(UChar_t) * kMaxRows);
 260     memset(&fClusters[0], 0x0, sizeof(AliTRDcluster*) * kMaxClustersLayer);
 261     memset(&fIndex[0], 0xffff, sizeof(UInt_t) * kMaxClustersLayer);
 262     return;
 263   }
 264
 265   // Make a copy
 266   AliTRDcluster *helpCL[kMaxClustersLayer];
 267   UInt_t helpInd[kMaxClustersLayer];
 268   nClStack = 0;
 269   // Defining iterators
 270   AliTRDcluster **fcliter = &fClusters[0], **hcliter = &helpCL[0]; UInt_t *finditer = &fIndex[0], *hinditer = &helpInd[0];
 271   AliTRDcluster *tmpcl = 0x0;
 272   for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
 273     if(!(tmpcl = *(fcliter++))){
 274         finditer++;
 275         continue;
 276     }
 277     *(hcliter++)  = tmpcl;
 278     *(hinditer++) = *finditer;
 279     tmpcl = 0x0;
 280     *(finditer++) = 0xffff;
 281     nClStack++;
 282   }
 283
 284   // do clusters arrangement
 285   fX = 0.;
 286   fN =  nClStack;
 287   nClStack = 0;
 288   // Reset Positions array
 289   memset(fPositions, 0, sizeof(UChar_t)*kMaxRows);
 290   AliTRDcluster **cliter = &helpCL[0]; // Declare iterator running over the whole array
 291   for(Int_t i = 0; i < fN; i++){
 292     // boundary check
 293     AliTRDcluster *cl = *(cliter++);
 294     UChar_t rowIndex = cl->GetPadRow();
 295     // Insert Leaf
 296     Int_t pos = FindYPosition(cl->GetY(), rowIndex, i);
 297     if(pos == -1){              // zbin is empty;
 298       Int_t upper = (rowIndex == fNRows - 1) ? nClStack : fPositions[rowIndex + 1];
 299       memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));
 300       memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
 301       fClusters[upper] = cl;
 302       fIndex[upper] = helpInd[i];
 303       // Move All pointer one position back
 304       for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
 305       nClStack++;
 306     } else {            // zbin not empty
 307       memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));
 308       memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
 309       fClusters[pos + 1] = cl;  //fIndex[i];
 310       fIndex[pos + 1] = helpInd[i];
 311       // Move All pointer one position back
 312       for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
 313       nClStack++;
 314     }
 315
 316     // calculate mean x
 317     fX += cl->GetX();
 318
 319     // Debug Streaming
 320     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 3){
 321       AliTRDcluster dcl(*cl);
 322       TTreeSRedirector &cstream = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
 323       cstream << "BuildIndices"
 324       << "Plane="    << fPlane
 325       << "Stack="    << fStack
 326       << "Sector="   << fSector
 327       << "Iter="     << iter
 328       << "C.="       << &dcl
 329       << "rowIndex=" << rowIndex
 330       << "\n";
 331     }
 332   }
 333
 334 //      AliInfo("Positions");
 335 //      for(int ir=0; ir<fNRows; ir++) printf("pos[%d] %d\n", ir, fPositions[ir]);
 336
 337   fX /= fN;
 338 }
 339
 340 //_____________________________________________________________________________
 341 Int_t AliTRDchamberTimeBin::Find(Float_t y) const
 342 {
 343   //
 344   // Returns index of the cluster nearest in Y
 345   //
 346
 347   if (fN <= 0) return 0;
 348
 349   if (y <= fClusters[0]->GetY()) return 0;
 350
 351   if (y >  fClusters[fN-1]->GetY()) return fN;
 352
 353
 354   Int_t b = 0;
 355   Int_t e = fN - 1;
 356   Int_t m = (b + e) / 2;
 357
 358   for ( ; b < e; m = (b + e) / 2) {
 359     if (y > fClusters[m]->GetY()) b = m + 1;
 360     else e = m;
 361   }
 362
 363   return m;
 364 }
 365
 366 //_____________________________________________________________________________
 367 Int_t AliTRDchamberTimeBin::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const
 368 {
 369 //
 370 // Tree search Algorithm to find the nearest left cluster for a given
 371 // y-position in a certain z-bin (in fact AVL-tree).
 372 // Making use of the fact that clusters are sorted in y-direction.
 373 //
 374 // Parameters:
 375 //   y : y position of the reference point in tracking coordinates
 376 //   z : z reference bin.
 377 //   nClusters :
 378 //
 379 // Output :
 380 // Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)
 381 //
 382
 383   Int_t start = fPositions[z];          // starting Position of the bin
 384   Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin
 385   Int_t end = upper - 1; // ending Position of the bin
 386   if(end < start) return -1; // Bin is empty
 387   Int_t middle = static_cast<Int_t>((start + end)/2);
 388   // 1st Part: climb down the tree: get the next cluster BEFORE ypos
 389   while(start + 1 < end){
 390     if(y >= fClusters[middle]->GetY()) start = middle;
 391     else end = middle;
 392     middle = static_cast<Int_t>((start + end)/2);
 393   }
 394   if(y > fClusters[end]->GetY()) return end;
 395   return start;
 396 }
 397
 398 //_____________________________________________________________________________
 399 Int_t AliTRDchamberTimeBin::FindNearestYCluster(Double_t y, UChar_t z) const
 400 {
 401 //
 402 // Tree search Algorithm to find the nearest cluster for a given
 403 // y-position in a certain z-bin (in fact AVL-tree).
 404 // Making use of the fact that clusters are sorted in y-direction.
 405 //
 406 // Parameters:
 407 //   y : y position of the reference point in tracking coordinates
 408 //   z : z reference bin.
 409 //
 410 // Output
 411 // Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)
 412 //
 413
 414   Int_t position = FindYPosition(y, z, fN);
 415   if(position == -1) return position; // bin empty
 416   // FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest
 417   // to the Reference y-position, so test both
 418   Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
 419   if((position + 1) < (upper)){
 420     if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
 421     else return position;
 422   }
 423   return position;
 424 }
 425
 426 //_____________________________________________________________________________
 427 Int_t AliTRDchamberTimeBin::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const
 428 {
 429 //
 430 // Finds the nearest cluster from a given point in a defined range.
 431 // Distance is determined in a 2D space by the 2-Norm.
 432 //
 433 // Parameters :
 434 //   y : y position of the reference point in tracking coordinates
 435 //   z : z reference bin.
 436 //   maxroady : maximum searching distance in y direction
 437 //   maxroadz : maximum searching distance in z direction
 438 //
 439 // Output
 440 // Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).
 441 // Cluster can be accessed with the operator[] or GetCluster(Int_t index)
 442 //
 443 // Detail description
 444 //
 445 // The following steps are perfomed:
 446 // 1. Get the expected z bins inside maxroadz.
 447 // 2. For each z bin find nearest y cluster.
 448 // 3. Select best candidate
 449 //
 450   Int_t   index   = -1;
 451   // initial minimal distance will be represented as ellipse: semi-major = z-direction
 452   // later 2-Norm will be used
 453 //      Float_t nExcentricity = TMath::Sqrt(maxroadz*maxroadz - maxroad*maxroad)/maxroadz;
 454   Float_t mindist = maxroadz;
 455
 456   // not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How
 457   // much neighbors depends on the Quotient maxroadz/fZLength
 458   UChar_t maxRows = 3;
 459   UChar_t zpos[kMaxRows];
 460   // Float_t mindist = TMath::Sqrt(maxroad*maxroad + maxroadz*maxroadz);
 461 //      UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);
 462   UChar_t myZbin = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
 463   if(z < fZ0) myZbin = fNRows - 1;
 464   if(z > fZ0 + fZLength) myZbin = 0;
 465   //printf("\n%f < %f < %f [%d]\n", fZ0, z, fZ0 + fZLength, myZbin);
 466   //for(int ic=0; ic<fN; ic++) printf("%d z = %f row %d\n", ic, fClusters[ic]->GetZ(), fClusters[ic]->GetPadRow());
 467
 468   UChar_t nNeighbors = 0;
 469   for(UChar_t i = 0; i < maxRows; i++){
 470     if((myZbin - 1 + i) < 0) continue;
 471     if((myZbin - 1 + i) > fNRows - 1) break;
 472     zpos[nNeighbors] = myZbin - 1 + i;
 473     nNeighbors++;
 474   }
 475   Float_t ycl = 0, zcl = 0;
 476   for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){ // Always test the neighbors too
 477     Int_t pos = FindNearestYCluster(y, zpos[neighbor]);
 478     if(pos == -1) continue;     // No cluster in bin
 479     AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);
 480     if(c->IsUsed()) continue;           // we are only interested in unused clusters
 481     ycl = c->GetY();
 482     // Too far away in y-direction (Prearrangement)
 483     if (TMath::Abs(ycl - y) > maxroady){
 484       //printf("y[%f] ycl[%f] roady[%f]\n", y, ycl, maxroady);
 485       continue;
 486     }
 487     zcl = c->GetZ();
 488     // Too far away in z-Direction
 489     // (Prearrangement since we have not so many bins to test)
 490     if (TMath::Abs(zcl - z) > maxroadz) continue;
 491
 492     Float_t dist;               // distance defined as 2-Norm
 493     // if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and
 494     // maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we
 495     // take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely
 496     // large enough to be usable as an indicator whether we have found a nearer cluster or not)
 497 //              if(mindist > 10000.){
 498 //                      Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));
 499 //                      mindist = maxroad/TMath::Sqrt(1 - nExcentricity*nExcentricity * (TMath::Cos(phi))*(TMath::Cos(phi)));
 500 //              }
 501     dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl));     // infinity Norm
 502 //              dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z));
 503     if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){
 504     //if((dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z))) < mindist){
 505       mindist = dist;
 506       index   = pos;
 507     }
 508   }
 509   // This is the Array Position in fIndex2D of the Nearest cluster: if a
 510   // cluster is called, then the function has to retrieve the Information
 511   // which is Stored in the Array called, the function
 512   return index;
 513 }
 514
 515 //_____________________________________________________________________________
 516 void AliTRDchamberTimeBin::BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta, Double_t phi)
 517 {
 518 // Helper function to calculate the area where to expect a cluster in THIS
 519 // layer.
 520 //
 521 // Parameters :
 522 //   cl    :
 523 //   cond  :
 524 //   Layer :
 525 //   theta :
 526 //   phi   :
 527 //
 528 // Detail description
 529 //
 530 // Helper function to calculate the area where to expect a cluster in THIS
 531 // layer. by using the information of a former cluster in another layer
 532 // and the angle in theta- and phi-direction between layer 0 and layer 3.
 533 // If the layer is zero, initial conditions are calculated. Otherwise a
 534 // linear interpolation is performed.
 535 //Begin_Html
 536 //<img src="gif/build_cond.gif">
 537 //End_Html
 538 //
 539
 540   if(!fReconstructor){
 541     AliError("Reconstructor not set.");
 542     return;
 543   }
 544
 545   if(Layer == 0){
 546     cond[0] = cl->GetY();                       // center: y-Direction
 547     cond[1] = cl->GetZ();                       // center: z-Direction
 548     cond[2] = fReconstructor->GetRecoParam()->GetMaxPhi()   * (cl->GetX() - GetX()) + 1.0;                      // deviation: y-Direction
 549     cond[3] = fReconstructor->GetRecoParam()->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0;                      // deviation: z-Direction
 550   } else {
 551     cond[0] = cl->GetY() + phi   * (GetX() - cl->GetX());
 552     cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
 553     cond[2] = fReconstructor->GetRecoParam()->GetRoad0y() + phi;
 554     cond[3] = fReconstructor->GetRecoParam()->GetRoad0z();
 555   }
 556 }
 557
 558 //_____________________________________________________________________________
 559 void AliTRDchamberTimeBin::GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize)
 560 {
 561 // Finds all clusters situated in this layer inside a rectangle  given by the center an ranges.
 562 //
 563 // Parameters :
 564 //   cond  :
 565 //   index :
 566 //   ncl :
 567 //   BufferSize   :
 568 //
 569 // Output :
 570 //
 571 // Detail description
 572 //
 573 // Function returs an array containing the indices in the stacklayer of
 574 // the clusters found an  the number of found clusters in the stacklayer
 575
 576   ncl = 0;
 577   memset(index, 0, BufferSize*sizeof(Int_t));
 578   if(fN == 0) return;
 579
 580   //Boundary checks
 581   Double_t zvals[2];
 582   if(((cond[1] - cond[3]) >= (fZ0 + fZLength)) || (cond[1] + cond[3]) <= fZ0) return; // We are outside of the chamvber
 583   zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
 584   zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength - 1.E-3;
 585
 586   UChar_t zhi = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
 587   UChar_t zlo = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
 588
 589 /*      AliInfo(Form("yc[%f] zc[%f] dy[%f] dz[%f]", cond[0], cond[1], cond[2], cond[3]));
 590   PrintClusters();
 591   AliInfo(Form("zlo[%f] zhi[%f]", zvals[0], zvals[1]));
 592   AliInfo(Form("zlo[%d] zhi[%d]", zlo, zhi));*/
 593
 594   Double_t ylo = cond[0] - cond[2],
 595            yhi = cond[0] + cond[2];
 596   //printf("CTB : ylo[%f] yhi[%f]\n", ylo, yhi);
 597   //Preordering in Direction z saves a lot of loops (boundary checked)
 598   for(UChar_t z = zlo; z <= zhi; z++){
 599     UInt_t upper = (z < fNRows-1) ? fPositions[z+1] : fN;
 600     //AliInfo(Form("z[%d] y [%d %d]", z, fPositions[z], upper));
 601     for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
 602       if(ncl == BufferSize){
 603         AliWarning("Buffer size riched. Some clusters may be lost.");
 604         return; //Buffer filled
 605       }
 606
 607       if(fClusters[y]->GetY() > yhi) break;                     // Abbortion conditions!!!
 608       if(fClusters[y]->GetY() < ylo) continue;  // Too small
 609       if(((Int_t)((fClusters[y]->GetZ())*1000) < (Int_t)(zvals[0]*1000)) || ((Int_t)((fClusters[y]->GetZ())*1000) > (Int_t)(zvals[1]*1000))){/*printf("exit z\n"); TODO*/ continue;}
 610       index[ncl] = y;
 611       ncl++;
 612     }
 613   }
 614   if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
 615 }
 616
 617 //_____________________________________________________________________________
 618 AliTRDcluster *AliTRDchamberTimeBin::GetNearestCluster(Double_t *cond)
 619 {
 620 // Function returning a pointer to the nearest cluster (nullpointer if not successfull).
 621 //
 622 // Parameters :
 623 //   cond  :
 624 //
 625 // Output :
 626 //   pointer to the nearest cluster (nullpointer if not successfull).
 627 //
 628 // Detail description
 629 //
 630 // returns a pointer to the nearest cluster (nullpointer if not
 631 // successfull) by the help of the method FindNearestCluster
 632
 633
 634   Double_t maxroad  = fReconstructor->GetRecoParam()->GetRoad2y();
 635   Double_t maxroadz = fReconstructor->GetRecoParam()->GetRoad2z();
 636
 637   Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
 638   AliTRDcluster *returnCluster = 0x0;
 639   if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
 640   return returnCluster;
 641 }
 642
 643 //_____________________________________________________________________________
 644 void AliTRDchamberTimeBin::Print(Option_t *) const
 645 {
 646 // Prints the position of each cluster in the stacklayer on the stdout
 647 //
 648   if(!fN) return;
 649   AliInfo(Form("Layer[%d] Stack[%d] Sector[%2d] nRows[%2d]", fPlane, fStack, fSector, fNRows));
 650   AliInfo(Form("Z0[%7.3f] Z1[%7.3f]", fZ0, fZ0+fZLength));
 651   AliTRDcluster* const* icl = &fClusters[0];
 652   for(Int_t jcl = 0; jcl < fN; jcl++, icl++){
 653     AliInfo(Form("%2d X[%7.3f] Y[%7.3f] Z[%7.3f] tb[%2d] col[%3d] row[%2d] used[%s]", jcl,  (*icl)->GetX(), (*icl)->GetY(), (*icl)->GetZ(), (*icl)->GetLocalTimeBin(), (*icl)->GetPadCol(), (*icl)->GetPadRow(),
 654     (*icl)->IsUsed() ? "y" : "n"));
 655   }
 656   Int_t irow = 0;
 657   for(UChar_t const* pos = &fPositions[0]; irow<fNRows; pos++, irow++){
 658     if((*pos) != 0xff) AliInfo(Form("r[%2d] pos[%d]", irow, (*pos)));
 659   }
 660 }