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