TRD/AliTRDchamberTimeBin.cxx

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