Line data Source code
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 48 : ClassImp(AliTRDchamberTimeBin)
46 :
47 : //_____________________________________________________________________________
48 : AliTRDchamberTimeBin::AliTRDchamberTimeBin(Int_t plane, Int_t stack, Int_t sector, Double_t z0, Double_t zLength)
49 10819 : :TObject()
50 10819 : ,fkReconstructor(NULL)
51 10819 : ,fPlane(plane)
52 10819 : ,fStack(stack)
53 10819 : ,fSector(sector)
54 10819 : ,fNRows(kMaxRows)
55 10819 : ,fN(0)
56 10819 : ,fX(0.)
57 10819 : ,fZ0(z0)
58 10819 : ,fZLength(zLength)
59 54095 : {
60 : //
61 : // Default constructor (Only provided to use AliTRDchamberTimeBin with arrays)
62 : //
63 10819 : SetBit(kT0, kFALSE);
64 10819 : SetBit(kOwner, kFALSE);
65 10819 : memset(fPositions, 1, kMaxRows*sizeof(UChar_t));
66 10819 : memset(fClusters, 0, kMaxClustersLayer*sizeof(AliTRDcluster*));
67 10819 : memset(fIndex, 1, kMaxClustersLayer*sizeof(UInt_t));
68 21638 : }
69 :
70 :
71 : //_____________________________________________________________________________
72 : AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer):
73 0 : TObject()
74 0 : ,fkReconstructor(layer.fkReconstructor)
75 0 : ,fPlane(layer.fPlane)
76 0 : ,fStack(layer.fStack)
77 0 : ,fSector(layer.fSector)
78 0 : ,fNRows(layer.fNRows)
79 0 : ,fN(layer.fN)
80 0 : ,fX(layer.fX)
81 0 : ,fZ0(layer.fZ0)
82 0 : ,fZLength(layer.fZLength)
83 0 : {
84 : // Copy Constructor
85 :
86 0 : SetBit(kT0, layer.IsT0());
87 0 : SetBit(kOwner, kFALSE);
88 0 : for(int i=0; i<kMaxRows; i++) fPositions[i] = layer.fPositions[i];
89 0 : memcpy(&fClusters[0], &layer.fClusters[0], kMaxClustersLayer*sizeof(AliTRDcluster*));
90 0 : memcpy(&fIndex[0], &layer.fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
91 :
92 :
93 : // BuildIndices();
94 0 : }
95 :
96 : //_____________________________________________________________________________
97 : AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDchamberTimeBin &layer)
98 : {
99 : // Assignment operator
100 :
101 0 : if (this != &layer) layer.Copy(*this);
102 0 : return *this;
103 : }
104 :
105 : //_____________________________________________________________________________
106 : void AliTRDchamberTimeBin::Clear(const Option_t *)
107 : {
108 : // Reset the Chamber Timebin
109 21638 : if(IsOwner())
110 0 : for(Int_t it = 0; it<kMaxClustersLayer; it++)
111 0 : delete fClusters[it];
112 10819 : memset(fClusters,0,kMaxClustersLayer*sizeof(fClusters[0]));
113 10819 : fN = 0;
114 10819 : }
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 0 : AliTRDchamberTimeBin &layer = (AliTRDchamberTimeBin &)o;
122 0 : layer.fkReconstructor = fkReconstructor;
123 0 : layer.fPlane = fPlane;
124 0 : layer.fStack = fStack;
125 0 : layer.fSector = fSector;
126 0 : layer.fNRows = fNRows;
127 0 : layer.fN = fN;
128 0 : layer.fX = fX;
129 0 : layer.fZ0 = fZ0;
130 0 : layer.fZLength = fZLength;
131 0 : layer.SetT0(IsT0());
132 0 : layer.SetBit(kOwner, kFALSE);
133 :
134 0 : for(int i=0; i<kMaxRows; i++) layer.fPositions[i] = fPositions[i];
135 0 : memcpy(&layer.fClusters[0], &fClusters[0], kMaxClustersLayer*sizeof(AliTRDcluster*));
136 0 : memcpy(&layer.fIndex[0], &fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
137 :
138 0 : TObject::Copy(layer); // copies everything into layer
139 :
140 : // layer.BuildIndices();
141 0 : }
142 :
143 : //_____________________________________________________________________________
144 : AliTRDchamberTimeBin::~AliTRDchamberTimeBin()
145 43276 : {
146 : // Destructor
147 10819 : if(IsOwner()){
148 0 : for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++) delete (*cit);
149 0 : }
150 21638 : }
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 0 : SetBit(kOwner, kTRUE);
161 0 : if(!copy) return;
162 0 : for(AliTRDcluster **cit = &fClusters[0]; (*cit); ++cit){
163 0 : (*cit) = new AliTRDcluster(*(*cit));
164 : }
165 0 : }
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 15594 : fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
177 7797 : fZLength = TMath::Abs(zLength);
178 7797 : }
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 36904 : if (fN == (Int_t) kMaxClustersLayer) {
194 : //AliWarning("Too many clusters !\n");
195 : return;
196 : }
197 :
198 18452 : if (fN == 0) {
199 7797 : fIndex[0] = index;
200 7797 : fClusters[fN++] = c;
201 7797 : return;
202 : }
203 :
204 10655 : Int_t i = Find(c->GetY());
205 10655 : memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
206 10655 : memmove(fIndex +i+1,fIndex +i,(fN-i)*sizeof(UInt_t));
207 10655 : fIndex[i] = index;
208 10655 : fClusters[i] = c;
209 10655 : fN++;
210 29107 : }
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 0 : fkReconstructor = rec;
219 0 : fPlane = AliTRDgeometry::GetLayer(det);
220 0 : fStack = AliTRDgeometry::GetStack(det);
221 0 : fSector = AliTRDgeometry::GetSector(det);
222 0 : AliTRDgeometry g;
223 0 : fNRows = g.GetPadPlane(fPlane, fStack)->GetNrows();
224 0 : fN = 0;
225 0 : for(AliTRDcluster **cit = &fClusters[0]; (*cit); cit++) fN++;
226 0 : BuildIndices();
227 0 : }
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 7797 : if(!fN) return;
244 :
245 : // Select clusters that belong to the Stack
246 : Int_t nClStack = 0; // Internal counter
247 52498 : for(Int_t i = 0; i < fN; i++){
248 36904 : if(fClusters[i]->IsUsed() || fClusters[i]->IsShared()){
249 0 : fClusters[i] = NULL;
250 0 : fIndex[i] = 0xffff;
251 18452 : } else nClStack++;
252 : }
253 7797 : 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 7797 : if(!nClStack){
257 0 : fN = 0;
258 0 : memset(&fPositions[0], 0, sizeof(UChar_t) * kMaxRows);
259 0 : memset(&fClusters[0], 0, sizeof(AliTRDcluster*) * kMaxClustersLayer);
260 0 : memset(&fIndex[0], 0, sizeof(UInt_t) * kMaxClustersLayer);
261 0 : return;
262 : }
263 :
264 : // Make a copy
265 7797 : AliTRDcluster *helpCL[kMaxClustersLayer];
266 7797 : UInt_t helpInd[kMaxClustersLayer];
267 : nClStack = 0;
268 : // Defining iterators
269 7797 : AliTRDcluster **fcliter = &fClusters[0], **hcliter = &helpCL[0]; UInt_t *finditer = &fIndex[0], *hinditer = &helpInd[0];
270 : AliTRDcluster *tmpcl = NULL;
271 52498 : for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
272 18452 : if(!(tmpcl = *(fcliter++))){
273 0 : finditer++;
274 0 : continue;
275 : }
276 18452 : *(hcliter++) = tmpcl;
277 18452 : *(hinditer++) = *finditer;
278 : tmpcl = NULL;
279 18452 : *(finditer++) = 0xffff;
280 18452 : nClStack++;
281 18452 : }
282 :
283 : // do clusters arrangement
284 7797 : fX = 0.;
285 7797 : fN = nClStack;
286 : nClStack = 0;
287 : // Reset Positions array
288 7797 : memset(fPositions, 0, sizeof(UChar_t)*kMaxRows);
289 : AliTRDcluster **cliter = &helpCL[0]; // Declare iterator running over the whole array
290 7797 : 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 52498 : for(Int_t i = 0; i < fN; i++){
293 : // boundary check
294 18452 : AliTRDcluster *cl = *(cliter++);
295 18452 : UChar_t rowIndex = cl->GetPadRow();
296 26249 : if(tb<0) tb=cl->GetLocalTimeBin();
297 : // Insert Leaf
298 18452 : Int_t pos = FindYPosition(cl->GetY(), rowIndex, nClStack);
299 18452 : if(pos == -2) continue; // zbin error;
300 18452 : else if(pos == -1) { // zbin is empty;
301 34185 : Int_t upper = (rowIndex == fNRows - 1) ? nClStack : fPositions[rowIndex + 1];
302 11395 : memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));
303 11395 : memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
304 11395 : fClusters[upper] = cl;
305 11395 : fIndex[upper] = helpInd[i];
306 : // Move All pointer one position back
307 172288 : for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
308 11395 : nClStack++;
309 11395 : } else { // zbin not empty
310 7057 : memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));
311 7057 : memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
312 7057 : fClusters[pos + 1] = cl; //fIndex[i];
313 7057 : fIndex[pos + 1] = helpInd[i];
314 : // Move All pointer one position back
315 86068 : for(UChar_t j = rowIndex + 1; j < fNRows; j++) fPositions[j]++;
316 7057 : nClStack++;
317 : }
318 :
319 : // calculate mean x
320 18452 : fX += cl->GetX();
321 :
322 : // Debug Streaming
323 18452 : if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) >= 3){
324 0 : AliTRDcluster dcl(*cl);
325 0 : TTreeSRedirector &cstream = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
326 0 : cstream << "BuildIndices"
327 0 : << "Plane=" << fPlane
328 0 : << "Stack=" << fStack
329 0 : << "Sector=" << fSector
330 0 : << "Iter=" << iter
331 0 : << "C.=" << &dcl
332 0 : << "rowIndex=" << rowIndex
333 0 : << "\n";
334 0 : }
335 36904 : }
336 :
337 : // AliInfo("Positions");
338 : // for(int ir=0; ir<fNRows; ir++) printf("pos[%d] %d\n", ir, fPositions[ir]);
339 7797 : if(nClStack < fN){
340 0 : AliWarning(Form("Found %d out of %d clusters outside in ChamberTimeBin[%02d_%d_%d|%2d]", fN-nClStack, fN, fSector, fStack, fPlane, tb));
341 0 : fN = nClStack;
342 0 : if(!fN){ // Nothing left in this time bin. Reset indexes
343 0 : memset(&fPositions[0], 0, sizeof(UChar_t) * kMaxRows);
344 0 : memset(&fClusters[0], 0, sizeof(AliTRDcluster*) * kMaxClustersLayer);
345 0 : memset(&fIndex[0], 0, sizeof(UInt_t) * kMaxClustersLayer);
346 0 : return;
347 : }
348 : }
349 15594 : if(fN) fX /= fN;
350 23391 : }
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 21310 : if (fN <= 0) return 0;
360 :
361 15205 : if (y <= fClusters[0]->GetY()) return 0;
362 :
363 7871 : if (y > fClusters[fN-1]->GetY()) return fN;
364 :
365 :
366 : Int_t b = 0;
367 4339 : Int_t e = fN - 1;
368 4339 : Int_t m = (b + e) / 2;
369 :
370 34356 : for ( ; b < e; m = (b + e) / 2) {
371 19554 : if (y > fClusters[m]->GetY()) b = m + 1;
372 : else e = m;
373 : }
374 :
375 : return m;
376 10655 : }
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 36904 : if(z>=fNRows){ // check pad row of cluster
396 0 : AliDebug(1, Form("Row[%2d] outside range [0 %2d] in %02d_%d_%d.", z, fNRows, fSector, fStack, fPlane));
397 0 : return -2;
398 : }
399 18452 : Int_t start = fPositions[z]; // starting Position of the bin
400 55356 : Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin
401 18452 : Int_t end = upper - 1; // ending Position of the bin
402 29847 : if(end < start) return -1; // Bin is empty
403 7057 : Int_t middle = static_cast<Int_t>((start + end)/2);
404 : // 1st Part: climb down the tree: get the next cluster BEFORE ypos
405 23942 : while(start + 1 < end){
406 9828 : if(y >= fClusters[middle]->GetY()) start = middle;
407 : else end = middle;
408 4914 : middle = static_cast<Int_t>((start + end)/2);
409 : }
410 14114 : if(y > fClusters[end]->GetY()) return end;
411 0 : return start;
412 18452 : }
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 0 : Int_t position = FindYPosition(y, z, fN);
431 0 : 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 0 : Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
435 0 : if((position + 1) < (upper)){
436 0 : if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
437 0 : else return position;
438 : }
439 0 : return position;
440 0 : }
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 0 : 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 0 : 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 0 : UChar_t myZbin = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
479 0 : if(z < fZ0) myZbin = fNRows - 1;
480 0 : 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 0 : for(UChar_t i = 0; i < maxRows; i++){
486 0 : if((myZbin - 1 + i) < 0) continue;
487 0 : if((myZbin - 1 + i) > fNRows - 1) break;
488 0 : zpos[nNeighbors] = myZbin - 1 + i;
489 0 : nNeighbors++;
490 0 : }
491 : Float_t ycl = 0, zcl = 0;
492 0 : for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){ // Always test the neighbors too
493 0 : Int_t pos = FindNearestYCluster(y, zpos[neighbor]);
494 0 : if(pos == -1) continue; // No cluster in bin
495 0 : AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);
496 0 : if(c->IsUsed()) continue; // we are only interested in unused clusters
497 0 : ycl = c->GetY();
498 : // Too far away in y-direction (Prearrangement)
499 0 : if (TMath::Abs(ycl - y) > maxroady){
500 : //printf("y[%f] ycl[%f] roady[%f]\n", y, ycl, maxroady);
501 0 : continue;
502 : }
503 0 : zcl = c->GetZ();
504 : // Too far away in z-Direction
505 : // (Prearrangement since we have not so many bins to test)
506 0 : 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 0 : 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 0 : 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 0 : }
524 0 : }
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 0 : return index;
529 0 : }
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 0 : if(!fkReconstructor){
557 0 : AliError("Reconstructor not set.");
558 0 : return;
559 : }
560 :
561 0 : if(Layer == 0){
562 0 : cond[0] = cl->GetY(); // center: y-Direction
563 0 : cond[1] = cl->GetZ(); // center: z-Direction
564 0 : cond[2] = fkReconstructor->GetRecoParam()->GetMaxPhi() * (cl->GetX() - GetX()) + 1.0; // deviation: y-Direction
565 0 : cond[3] = fkReconstructor->GetRecoParam()->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0; // deviation: z-Direction
566 0 : } else {
567 0 : cond[0] = cl->GetY() + phi * (GetX() - cl->GetX());
568 0 : cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
569 0 : cond[2] = fkReconstructor->GetRecoParam()->GetRoad0y() + phi;
570 0 : cond[3] = fkReconstructor->GetRecoParam()->GetRoad0z();
571 : }
572 0 : }
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 11790 : ncl = 0;
593 5895 : memset(index, 0, BufferSize*sizeof(Int_t));
594 5895 : if(fN == 0) return;
595 :
596 : //Boundary checks
597 : Double_t zvals[2];
598 11790 : if(((cond[1] - cond[3]) >= (fZ0 + fZLength)) || (cond[1] + cond[3]) <= fZ0) return; // We are outside of the chamvber
599 17685 : zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
600 17685 : zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength - 1.E-3;
601 :
602 5895 : UChar_t zhi = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
603 5895 : 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 5895 : Double_t ylo = cond[0] - cond[2],
611 5895 : 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 99835 : for(UChar_t z = zlo; z <= zhi; z++){
615 123225 : UInt_t upper = (z < fNRows-1) ? fPositions[z+1] : fN;
616 : //AliInfo(Form("z[%d] y [%d %d]", z, fPositions[z], upper));
617 138135 : for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
618 8496 : if(ncl == BufferSize){
619 0 : AliDebug(1, Form("Buffer size [%d] riched. Some clusters may be lost.", BufferSize));
620 0 : return; //Buffer filled
621 : }
622 :
623 9190 : if(fClusters[y]->GetY() > yhi) break; // Abbortion conditions!!!
624 7802 : if(fClusters[y]->GetY() < ylo) continue; // Too small
625 11062 : 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 5487 : index[ncl] = y;
627 5487 : ncl++;
628 5487 : }
629 41075 : }
630 5895 : if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
631 11790 : }
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 0 : Double_t maxroad = fkReconstructor->GetRecoParam()->GetRoad2y();
651 0 : Double_t maxroadz = fkReconstructor->GetRecoParam()->GetRoad2z();
652 :
653 0 : Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
654 : AliTRDcluster *returnCluster = NULL;
655 0 : if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
656 0 : 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 0 : if(!fN) return;
665 0 : AliInfo(Form("Layer[%d] Stack[%d] Sector[%2d] nRows[%2d]", fPlane, fStack, fSector, fNRows));
666 0 : AliInfo(Form("Z0[%7.3f] Z1[%7.3f]", fZ0, fZ0+fZLength));
667 0 : AliTRDcluster* const* icl = &fClusters[0];
668 0 : for(Int_t jcl = 0; jcl < fN; jcl++, icl++){
669 0 : 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 0 : for(UChar_t const* pos = &fPositions[0]; irow<fNRows; pos++, irow++){
674 0 : if((*pos) != 0xff) AliInfo(Form("r[%2d] pos[%d]", irow, (*pos)));
675 : }
676 0 : }
|
