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