1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDchamberTimeBin.cxx 23313 2008-01-11 14:56:43Z cblume $ */
18 ///////////////////////////////////////////////////////////////////////////////
20 // Organization of clusters at the level of 1 TRD chamber. //
21 // The data structure is used for tracking at the stack level. //
23 // Functionalities: //
24 // 1. cluster organization and sorting //
25 // 2. fast data navigation //
28 // Alex Bercuci <A.Bercuci@gsi.de> //
29 // Markus Fasel <M.Fasel@gsi.de> //
31 ///////////////////////////////////////////////////////////////////////////////
36 #include <TStopwatch.h>
37 #include <TTreeStream.h>
41 #include "AliTRDgeometry.h"
42 #include "AliTRDpadPlane.h"
43 #include "AliTRDchamberTimeBin.h"
44 #include "AliTRDrecoParam.h"
45 #include "AliTRDReconstructor.h"
46 #include "AliTRDtrackerV1.h"
49 ClassImp(AliTRDchamberTimeBin)
51 //_____________________________________________________________________________
52 AliTRDchamberTimeBin::AliTRDchamberTimeBin(Int_t plane, Int_t stack, Int_t sector, Double_t z0, Double_t zLength)
65 // Default constructor (Only provided to use AliTRDchamberTimeBin with arrays)
68 SetBit(kOwner, kFALSE);
69 memset(fPositions, 1, kMaxRows*sizeof(UChar_t));
70 memset(fClusters, 0, kMaxClustersLayer*sizeof(AliTRDcluster*));
71 memset(fIndex, 1, kMaxClustersLayer*sizeof(UInt_t));
75 //_____________________________________________________________________________
76 AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer):
78 ,fReconstructor(layer.fReconstructor)
81 ,fSector(layer.fSector)
86 ,fZLength(layer.fZLength)
90 SetBit(kT0, layer.IsT0());
91 SetBit(kOwner, kFALSE);
92 for(int i=0; i<kMaxRows; i++) fPositions[i] = layer.fPositions[i];
93 memcpy(&fClusters[0], &layer.fClusters[0], kMaxClustersLayer*sizeof(AliTRDcluster*));
94 memcpy(&fIndex[0], &layer.fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
100 //_____________________________________________________________________________
101 AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDchamberTimeBin &layer)
103 // Assignment operator
105 if (this != &layer) layer.Copy(*this);
109 //_____________________________________________________________________________
110 void AliTRDchamberTimeBin::Clear(const Option_t *)
112 for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++){
113 if(IsOwner()) delete (*cit);
119 //_____________________________________________________________________________
120 void AliTRDchamberTimeBin::Copy(TObject &o) const
122 // Copy method. Performs a deep copy of all data from this object to object o.
124 AliTRDchamberTimeBin &layer = (AliTRDchamberTimeBin &)o;
125 layer.fReconstructor = fReconstructor;
126 layer.fPlane = fPlane;
127 layer.fStack = fStack;
128 layer.fSector = fSector;
129 layer.fNRows = fNRows;
133 layer.fZLength = fZLength;
135 layer.SetBit(kOwner, kFALSE);
137 for(int i=0; i<kMaxRows; i++) layer.fPositions[i] = fPositions[i];
138 memcpy(&layer.fClusters[0], &fClusters[0], kMaxClustersLayer*sizeof(AliTRDcluster*));
139 memcpy(&layer.fIndex[0], &fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
141 TObject::Copy(layer); // copies everything into layer
143 // layer.BuildIndices();
146 //_____________________________________________________________________________
147 AliTRDchamberTimeBin::~AliTRDchamberTimeBin()
151 for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++) delete (*cit);
155 //_____________________________________________________________________________
156 void AliTRDchamberTimeBin::SetOwner(Bool_t copy)
158 // Sets the ownership of the clusters to this
159 // If option "copy" is kTRUE [default] the clusters
160 // are also copied otherwise only the ownership bit
163 SetBit(kOwner, kTRUE);
165 for(AliTRDcluster **cit = &fClusters[0]; (*cit); ++cit){
166 (*cit) = new AliTRDcluster(*(*cit));
170 //_____________________________________________________________________________
171 void AliTRDchamberTimeBin::SetRange(Float_t z0, Float_t zLength)
173 // Sets the range in z-direction
176 // z0 : starting position of layer in the z direction
177 // zLength : length of layer in the z direction
179 fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
180 fZLength = TMath::Abs(zLength);
183 //_____________________________________________________________________________
184 void AliTRDchamberTimeBin::InsertCluster(AliTRDcluster *c, UInt_t index)
187 // Insert cluster in cluster array.
188 // Clusters are sorted according to Y coordinate.
191 //if (fTimeBinIndex < 0) {
192 //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
196 if (fN == (Int_t) kMaxClustersLayer) {
197 //AliWarning("Too many clusters !\n");
207 Int_t i = Find(c->GetY());
208 memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
209 memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
215 //___________________________________________________
216 void AliTRDchamberTimeBin::Bootstrap(const AliTRDReconstructor *rec, Int_t det)
218 // Reinitialize all data members from the clusters array
219 // It has to be used after reading from disk
221 fReconstructor = rec;
222 fPlane = AliTRDgeometry::GetLayer(det);
223 fStack = AliTRDgeometry::GetStack(det);
224 fSector = AliTRDgeometry::GetSector(det);
226 fNRows = g.GetPadPlane(fPlane, fStack)->GetNrows();
228 for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++) fN++;
232 //_____________________________________________________________________________
233 void AliTRDchamberTimeBin::BuildIndices(Int_t iter)
235 // Rearrangement of the clusters belonging to the propagation layer for the stack.
237 // Detailed description
239 // The array indices of all clusters in one PropagationLayer are stored in
240 // array. The array is divided into several bins.
241 // The clusters are sorted in increasing order of their y coordinate.
243 // Sorting algorithm: TreeSearch
248 // Select clusters that belong to the Stack
249 Int_t nClStack = 0; // Internal counter
250 for(Int_t i = 0; i < fN; i++){
251 if(fClusters[i]->IsUsed()){
256 if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d. Truncating.", nClStack, kMaxClustersLayer));
258 // Nothing in this time bin. Reset indexes
261 memset(&fPositions[0], 0xff, sizeof(UChar_t) * kMaxRows);
262 memset(&fClusters[0], 0x0, sizeof(AliTRDcluster*) * kMaxClustersLayer);
263 memset(&fIndex[0], 0xffff, sizeof(UInt_t) * kMaxClustersLayer);
268 AliTRDcluster *helpCL[kMaxClustersLayer];
269 UInt_t helpInd[kMaxClustersLayer];
271 // Defining iterators
272 AliTRDcluster **fcliter = &fClusters[0], **hcliter = &helpCL[0]; UInt_t *finditer = &fIndex[0], *hinditer = &helpInd[0];
273 AliTRDcluster *tmpcl = 0x0;
274 for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
275 if(!(tmpcl = *(fcliter++))){
279 *(hcliter++) = tmpcl;
280 *(hinditer++) = *finditer;
282 *(finditer++) = 0xffff;
286 // do clusters arrangement
290 // Reset Positions array
291 memset(fPositions, 0, sizeof(UChar_t)*kMaxRows);
292 AliTRDcluster **cliter = &helpCL[0]; // Declare iterator running over the whole array
293 for(Int_t i = 0; i < fN; i++){
295 AliTRDcluster *cl = *(cliter++);
296 UChar_t rowIndex = cl->GetPadRow();
298 Int_t pos = FindYPosition(cl->GetY(), rowIndex, i);
299 if(pos == -1){ // zbin is empty;
300 Int_t upper = (rowIndex == fNRows - 1) ? nClStack : fPositions[rowIndex + 1];
301 memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));
302 memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
303 fClusters[upper] = cl;
304 fIndex[upper] = helpInd[i];
305 // Move All pointer one position back
306 for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
308 } else { // zbin not empty
309 memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));
310 memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
311 fClusters[pos + 1] = cl; //fIndex[i];
312 fIndex[pos + 1] = helpInd[i];
313 // Move All pointer one position back
314 for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
322 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 3){
323 TTreeSRedirector &cstream = *AliTRDtrackerV1::DebugStreamer();
324 cstream << "BuildIndices"
325 << "Plane=" << fPlane
326 << "Stack=" << fStack
327 << "Sector=" << fSector
330 << "rowIndex=" << rowIndex
335 // AliInfo("Positions");
336 // for(int ir=0; ir<fNRows; ir++) printf("pos[%d] %d\n", ir, fPositions[ir]);
341 //_____________________________________________________________________________
342 Int_t AliTRDchamberTimeBin::Find(Float_t y) const
345 // Returns index of the cluster nearest in Y
348 if (fN <= 0) return 0;
350 if (y <= fClusters[0]->GetY()) return 0;
352 if (y > fClusters[fN-1]->GetY()) return fN;
357 Int_t m = (b + e) / 2;
359 for ( ; b < e; m = (b + e) / 2) {
360 if (y > fClusters[m]->GetY()) b = m + 1;
367 //_____________________________________________________________________________
368 Int_t AliTRDchamberTimeBin::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const
371 // Tree search Algorithm to find the nearest left cluster for a given
372 // y-position in a certain z-bin (in fact AVL-tree).
373 // Making use of the fact that clusters are sorted in y-direction.
376 // y : y position of the reference point in tracking coordinates
377 // z : z reference bin.
381 // Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)
384 Int_t start = fPositions[z]; // starting Position of the bin
385 Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin
386 Int_t end = upper - 1; // ending Position of the bin
387 if(end < start) return -1; // Bin is empty
388 Int_t middle = static_cast<Int_t>((start + end)/2);
389 // 1st Part: climb down the tree: get the next cluster BEFORE ypos
390 while(start + 1 < end){
391 if(y >= fClusters[middle]->GetY()) start = middle;
393 middle = static_cast<Int_t>((start + end)/2);
395 if(y > fClusters[end]->GetY()) return end;
399 //_____________________________________________________________________________
400 Int_t AliTRDchamberTimeBin::FindNearestYCluster(Double_t y, UChar_t z) const
403 // Tree search Algorithm to find the nearest cluster for a given
404 // y-position in a certain z-bin (in fact AVL-tree).
405 // Making use of the fact that clusters are sorted in y-direction.
408 // y : y position of the reference point in tracking coordinates
409 // z : z reference bin.
412 // Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)
415 Int_t position = FindYPosition(y, z, fN);
416 if(position == -1) return position; // bin empty
417 // FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest
418 // to the Reference y-position, so test both
419 Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
420 if((position + 1) < (upper)){
421 if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
422 else return position;
427 //_____________________________________________________________________________
428 Int_t AliTRDchamberTimeBin::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const
431 // Finds the nearest cluster from a given point in a defined range.
432 // Distance is determined in a 2D space by the 2-Norm.
435 // y : y position of the reference point in tracking coordinates
436 // z : z reference bin.
437 // maxroady : maximum searching distance in y direction
438 // maxroadz : maximum searching distance in z direction
441 // Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).
442 // Cluster can be accessed with the operator[] or GetCluster(Int_t index)
444 // Detail description
446 // The following steps are perfomed:
447 // 1. Get the expected z bins inside maxroadz.
448 // 2. For each z bin find nearest y cluster.
449 // 3. Select best candidate
452 // initial minimal distance will be represented as ellipse: semi-major = z-direction
453 // later 2-Norm will be used
454 // Float_t nExcentricity = TMath::Sqrt(maxroadz*maxroadz - maxroad*maxroad)/maxroadz;
455 Float_t mindist = maxroadz;
457 // not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How
458 // much neighbors depends on the Quotient maxroadz/fZLength
460 UChar_t zpos[kMaxRows];
461 // Float_t mindist = TMath::Sqrt(maxroad*maxroad + maxroadz*maxroadz);
462 // UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);
463 UChar_t myZbin = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
464 if(z < fZ0) myZbin = fNRows - 1;
465 if(z > fZ0 + fZLength) myZbin = 0;
466 //printf("\n%f < %f < %f [%d]\n", fZ0, z, fZ0 + fZLength, myZbin);
467 //for(int ic=0; ic<fN; ic++) printf("%d z = %f row %d\n", ic, fClusters[ic]->GetZ(), fClusters[ic]->GetPadRow());
469 UChar_t nNeighbors = 0;
470 for(UChar_t i = 0; i < maxRows; i++){
471 if((myZbin - 1 + i) < 0) continue;
472 if((myZbin - 1 + i) > fNRows - 1) break;
473 zpos[nNeighbors] = myZbin - 1 + i;
476 Float_t ycl = 0, zcl = 0;
477 for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){ // Always test the neighbors too
478 Int_t pos = FindNearestYCluster(y, zpos[neighbor]);
479 if(pos == -1) continue; // No cluster in bin
480 AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);
481 if(c->IsUsed()) continue; // we are only interested in unused clusters
483 // Too far away in y-direction (Prearrangement)
484 if (TMath::Abs(ycl - y) > maxroady){
485 //printf("y[%f] ycl[%f] roady[%f]\n", y, ycl, maxroady);
489 // Too far away in z-Direction
490 // (Prearrangement since we have not so many bins to test)
491 if (TMath::Abs(zcl - z) > maxroadz) continue;
493 Float_t dist; // distance defined as 2-Norm
494 // if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and
495 // maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we
496 // take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely
497 // large enough to be usable as an indicator whether we have found a nearer cluster or not)
498 // if(mindist > 10000.){
499 // Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));
500 // mindist = maxroad/TMath::Sqrt(1 - nExcentricity*nExcentricity * (TMath::Cos(phi))*(TMath::Cos(phi)));
502 dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl)); // infinity Norm
503 // dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z));
504 if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){
505 //if((dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z))) < mindist){
510 // This is the Array Position in fIndex2D of the Nearest cluster: if a
511 // cluster is called, then the function has to retrieve the Information
512 // which is Stored in the Array called, the function
516 //_____________________________________________________________________________
517 void AliTRDchamberTimeBin::BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta, Double_t phi)
519 // Helper function to calculate the area where to expect a cluster in THIS
529 // Detail description
531 // Helper function to calculate the area where to expect a cluster in THIS
532 // layer. by using the information of a former cluster in another layer
533 // and the angle in theta- and phi-direction between layer 0 and layer 3.
534 // If the layer is zero, initial conditions are calculated. Otherwise a
535 // linear interpolation is performed.
537 //<img src="gif/build_cond.gif">
542 AliError("Reconstructor not set.");
547 cond[0] = cl->GetY(); // center: y-Direction
548 cond[1] = cl->GetZ(); // center: z-Direction
549 cond[2] = fReconstructor->GetRecoParam()->GetMaxPhi() * (cl->GetX() - GetX()) + 1.0; // deviation: y-Direction
550 cond[3] = fReconstructor->GetRecoParam()->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0; // deviation: z-Direction
552 cond[0] = cl->GetY() + phi * (GetX() - cl->GetX());
553 cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
554 cond[2] = fReconstructor->GetRecoParam()->GetRoad0y() + phi;
555 cond[3] = fReconstructor->GetRecoParam()->GetRoad0z();
559 //_____________________________________________________________________________
560 void AliTRDchamberTimeBin::GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize)
562 // Finds all clusters situated in this layer inside a rectangle given by the center an ranges.
572 // Detail description
574 // Function returs an array containing the indices in the stacklayer of
575 // the clusters found an the number of found clusters in the stacklayer
578 memset(index, 0, BufferSize*sizeof(Int_t));
583 if(((cond[1] - cond[3]) >= (fZ0 + fZLength)) || (cond[1] + cond[3]) <= fZ0) return; // We are outside of the chamvber
584 zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
585 zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength - 1.E-3;
587 UChar_t zhi = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
588 UChar_t zlo = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
590 /* AliInfo(Form("yc[%f] zc[%f] dy[%f] dz[%f]", cond[0], cond[1], cond[2], cond[3]));
592 AliInfo(Form("zlo[%f] zhi[%f]", zvals[0], zvals[1]));
593 AliInfo(Form("zlo[%d] zhi[%d]", zlo, zhi));*/
595 //Preordering in Direction z saves a lot of loops (boundary checked)
596 for(UChar_t z = zlo; z <= zhi; z++){
597 UInt_t upper = (z < fNRows-1) ? fPositions[z+1] : fN;
598 //AliInfo(Form("z[%d] y [%d %d]", z, fPositions[z], upper));
599 for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
600 if(ncl == BufferSize){
601 AliWarning("Buffer size riched. Some clusters may be lost.");
602 return; //Buffer filled
605 if(fClusters[y]->GetY() > (cond[0] + cond[2])) break; // Abbortion conditions!!!
606 if(fClusters[y]->GetY() < (cond[0] - cond[2])) continue; // Too small
607 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;}
612 if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
615 //_____________________________________________________________________________
616 AliTRDcluster *AliTRDchamberTimeBin::GetNearestCluster(Double_t *cond)
618 // Function returning a pointer to the nearest cluster (nullpointer if not successfull).
624 // pointer to the nearest cluster (nullpointer if not successfull).
626 // Detail description
628 // returns a pointer to the nearest cluster (nullpointer if not
629 // successfull) by the help of the method FindNearestCluster
632 Double_t maxroad = fReconstructor->GetRecoParam()->GetRoad2y();
633 Double_t maxroadz = fReconstructor->GetRecoParam()->GetRoad2z();
635 Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
636 AliTRDcluster *returnCluster = 0x0;
637 if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
638 return returnCluster;
641 //_____________________________________________________________________________
642 void AliTRDchamberTimeBin::Print(Option_t *) const
644 // Prints the position of each cluster in the stacklayer on the stdout
647 AliInfo(Form("Layer[%d] Stack[%d] Sector[%2d] nRows[%2d]", fPlane, fStack, fSector, fNRows));
648 AliInfo(Form("Z0[%7.3f] Z1[%7.3f]", fZ0, fZ0+fZLength));
649 AliTRDcluster* const* icl = &fClusters[0];
650 for(Int_t jcl = 0; jcl < fN; jcl++, icl++){
651 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(),
652 (*icl)->IsUsed() ? "y" : "n"));
655 for(UChar_t const* pos = &fPositions[0]; irow<fNRows; pos++, irow++){
656 if((*pos) != 0xff) AliInfo(Form("r[%2d] pos[%d]", irow, (*pos)));