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