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: AliTRDtrackingChamber.cxx 23810 2008-02-08 09:00:27Z hristov $ */
17 :
18 : ////////////////////////////////////////////////////////////////////////////
19 : // //
20 : // Tracking in one chamber //
21 : // //
22 : // Authors: //
23 : // Alex Bercuci <A.Bercuci@gsi.de> //
24 : // Markus Fasel <M.Fasel@gsi.de> //
25 : // //
26 : ////////////////////////////////////////////////////////////////////////////
27 :
28 : #include "AliTRDtrackingChamber.h"
29 :
30 : #include "TMath.h"
31 : #include "TMatrixTBase.h"
32 : #include <TTreeStream.h>
33 :
34 : #include "AliTRDReconstructor.h"
35 : #include "AliTRDrecoParam.h"
36 : #include "AliTRDtrackerV1.h"
37 : #include "AliTRDgeometry.h"
38 : #include "AliTRDpadPlane.h"
39 : #include "AliTRDcalibDB.h"
40 : #include "AliTRDCommonParam.h"
41 : #include "AliTRDCalDet.h"
42 : #include "AliTRDCalROC.h"
43 :
44 48 : ClassImp(AliTRDtrackingChamber)
45 :
46 : //_______________________________________________________
47 21987 : AliTRDtrackingChamber::AliTRDtrackingChamber()
48 349 : :TObject()
49 349 : ,fDetector(-1)
50 349 : ,fX0(0.)
51 : // ,fExB(0.)
52 : // ,fVD(0.)
53 : // ,fT0(0.)
54 : // ,fS2PRF(0.)
55 : // ,fDiffL(0.)
56 : // ,fDiffT(0.)
57 2094 : {}
58 :
59 : //_______________________________________________________
60 : void AliTRDtrackingChamber::Clear(const Option_t *opt)
61 : {
62 22685 : for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++) fTB[itb].Clear(opt);
63 349 : }
64 :
65 : //_______________________________________________________
66 : Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *const geo, Bool_t hlt)
67 : {
68 : // Init chamber and all time bins (AliTRDchamberTimeBin)
69 : // Calculates radial position of the chamber based on
70 : // radial positions of the time bins (calibration/alignment aware)
71 : //
72 698 : if(fDetector < 0 || fDetector >= AliTRDgeometry::kNdet){
73 0 : AliWarning(Form("Detector index not set correctly to %d", fDetector));
74 0 : return kFALSE;
75 : }
76 :
77 349 : Int_t stack = AliTRDgeometry::GetStack(fDetector);
78 349 : Int_t layer = AliTRDgeometry::GetLayer(fDetector);
79 349 : AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
80 349 : Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
81 349 : Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;
82 349 : Int_t nrows = pp->GetNrows();
83 :
84 349 : Int_t index[50], jtb = 0;
85 22336 : for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
86 10819 : if(!fTB[itb]) continue;
87 7797 : fTB[itb].SetRange(z0, zl);
88 7797 : fTB[itb].SetNRows(nrows);
89 7797 : fTB[itb].SetPlane(layer);
90 7797 : fTB[itb].SetStack(stack);
91 7797 : fTB[itb].SetSector(AliTRDgeometry::GetSector(fDetector));
92 7797 : fTB[itb].BuildIndices();
93 7797 : index[jtb++] = itb;
94 7797 : }
95 351 : if(jtb<2) return kFALSE;
96 :
97 347 : AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
98 : Float_t t0;
99 347 : if(!hlt){
100 347 : t0 = calib->GetT0Average(fDetector);
101 347 : }else{
102 0 : t0 = calib->GetT0Det()->GetValue(fDetector);
103 : }
104 : // fVD = calib->GetVdriftAverage(fDetector);
105 : // fS2PRF = calib->GetPRFROC(fDetector)->GetMean(); fS2PRF *= fS2PRF;
106 : // fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
107 : // AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL, fDiffT, fVD);
108 :
109 : // ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
110 : //fTB[Int_t(t0)].SetT0();
111 347 : Double_t x0 = fTB[index[0]].GetX();
112 347 : Double_t x1 = fTB[index[1]].GetX();
113 347 : Double_t dx = (x0 - x1)/(index[1] - index[0]);
114 347 : fX0 = x0 + dx*(index[0] - t0);
115 : return kTRUE;
116 698 : }
117 :
118 : //_______________________________________________________
119 : Int_t AliTRDtrackingChamber::GetNClusters() const
120 : {
121 : // Basic loop method
122 : // Returns number of clusters in chamber
123 : //
124 : Int_t n = 0;
125 17160 : for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
126 8184 : n += Int_t(fTB[itb]);
127 : }
128 264 : return n;
129 : }
130 :
131 : //_______________________________________________________
132 : void AliTRDtrackingChamber::Bootstrap(const AliTRDReconstructor *rec)
133 : {
134 : // Basic loop method
135 : // Bootstrap each time bin
136 : //
137 0 : AliTRDchamberTimeBin *jtb = &fTB[0];
138 0 : for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){
139 0 : (*jtb).Bootstrap(rec, fDetector);
140 : }
141 0 : }
142 :
143 : //_______________________________________________________
144 : void AliTRDtrackingChamber::SetOwner()
145 : {
146 : // Basic loop method
147 : // Set ownership in time bins
148 : //
149 0 : AliTRDchamberTimeBin *jtb = &fTB[0];
150 0 : for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++, ++jtb){
151 0 : if(!(Int_t(*jtb))) continue;
152 0 : (*jtb).SetOwner();
153 0 : }
154 0 : }
155 :
156 : //_______________________________________________________
157 : Double_t AliTRDtrackingChamber::GetQuality()
158 : {
159 : //
160 : // Calculate chamber quality for seeding.
161 : //
162 : //
163 : // Parameters :
164 : // layers : Array of propagation layers for this plane.
165 : //
166 : // Output :
167 : // plane quality factor for seeding
168 : //
169 : // Detailed description
170 : //
171 : // The quality of the plane for seeding is higher if:
172 : // 1. the average timebin population is closer to an integer number
173 : // 2. the distribution of clusters/timebin is closer to a uniform distribution.
174 : // - the slope of the first derivative of a parabolic fit is small or
175 : // - the slope of a linear fit is small
176 : //
177 :
178 : Int_t ncl = 0;
179 : Int_t nused = 0;
180 : Int_t nClLayer;
181 0 : for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
182 0 : if(!(nClLayer = fTB[itb].GetNClusters())) continue;
183 0 : ncl += nClLayer;
184 0 : for(Int_t incl = 0; incl < nClLayer; incl++){
185 0 : if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
186 : }
187 0 : }
188 :
189 : // calculate the deviation of the mean number of clusters from the
190 : // closest integer values
191 0 : Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
192 0 : Int_t ncli = Int_t(nclMed);
193 0 : Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
194 0 : nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
195 0 : return TMath::Exp(-5.*TMath::Abs(nclDev));
196 :
197 : // // get slope of the derivative
198 : // if(!fitter.Eval()) return quality;
199 : // fitter.PrintResults(3);
200 : // Double_t a = fitter.GetParameter(1);
201 : //
202 : // printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);
203 : // return quality*TMath::Exp(-a);
204 :
205 : }
206 :
207 :
208 : //_______________________________________________________
209 : Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry * const geo, const AliTRDReconstructor *rec)
210 : {
211 : //
212 : // Creates a seeding layer
213 : //
214 :
215 : // constants
216 : const Int_t kMaxRows = 16;
217 : const Int_t kMaxCols = 144;
218 : const Int_t kMaxPads = 2304;
219 0 : Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();
220 0 : Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();
221 :
222 : // Get the geometrical data of the chamber
223 0 : Int_t layer = geo->GetLayer(fDetector);
224 0 : Int_t stack = geo->GetStack(fDetector);
225 0 : Int_t sector= geo->GetSector(fDetector);
226 0 : AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
227 0 : Int_t nCols = pp->GetNcols();
228 0 : Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
229 0 : Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
230 0 : Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
231 0 : Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
232 : Float_t z0 = -1., zl = -1.;
233 0 : Int_t nRows = pp->GetNrows();
234 0 : Float_t binlength = (ymax - ymin)/nCols;
235 : //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
236 :
237 : // Fill the histogram
238 : Int_t nClusters;
239 0 : Int_t *histogram[kMaxRows]; // 2D-Histogram
240 0 : Int_t hvals[kMaxPads + 1]; memset(hvals, 0, sizeof(Int_t)*kMaxPads); // one entry in addition for termination flag
241 0 : Float_t *sigmas[kMaxRows];
242 0 : Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads);
243 : AliTRDcluster *c = NULL;
244 0 : for(Int_t irs = 0; irs < kMaxRows; irs++){
245 0 : histogram[irs] = &hvals[irs*kMaxCols];
246 0 : sigmas[irs] = &svals[irs*kMaxCols];
247 : }
248 0 : for(Int_t iTime = timeBinMin; iTime < AliTRDseedV1::kNtb-timeBinMax; iTime++){
249 0 : if(!(nClusters = fTB[iTime].GetNClusters())) continue;
250 0 : z0 = fTB[iTime].GetZ0();
251 0 : zl = fTB[iTime].GetDZ0();
252 0 : for(Int_t incl = 0; incl < nClusters; incl++){
253 0 : c = fTB[iTime].GetCluster(incl);
254 0 : histogram[c->GetPadRow()][c->GetPadCol()]++;
255 0 : sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
256 : }
257 0 : }
258 :
259 : // Now I have everything in the histogram, do the selection
260 : //Int_t nPads = nCols * nRows;
261 : // This is what we are interested in: The center of gravity of the best candidates
262 0 : Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
263 0 : Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
264 0 : Float_t *cogy[kMaxRows];
265 0 : Float_t *cogz[kMaxRows];
266 :
267 : // Lookup-Table storing coordinates according to the bins
268 0 : Float_t yLengths[kMaxCols]; memset(yLengths, 0, kMaxCols*sizeof(Float_t));
269 0 : Float_t zLengths[kMaxRows]; memset(zLengths, 0, kMaxRows*sizeof(Float_t));
270 0 : for(Int_t icnt = 0; icnt < nCols; icnt++){
271 0 : yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
272 : }
273 0 : for(Int_t icnt = 0; icnt < nRows; icnt++){
274 0 : zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
275 : }
276 :
277 : // A bitfield is used to mask the pads as usable
278 0 : Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
279 0 : for(UChar_t icount = 0; icount < nRows; icount++){
280 0 : cogy[icount] = &cogyvals[icount*kMaxCols];
281 0 : cogz[icount] = &cogzvals[icount*kMaxCols];
282 : }
283 : // In this array the array position of the best candidates will be stored
284 0 : Int_t cand[AliTRDtrackerV1::kMaxTracksStack];
285 0 : Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
286 :
287 : // helper variables
288 0 : Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
289 : Int_t nCandidates = 0;
290 : Float_t norm, cogv;
291 : // histogram filled -> Select best bins
292 0 : Int_t nPads = nCols * nRows;
293 : // take out all the bins which have less than 3 entries (faster sorting)
294 0 : Int_t content[kMaxPads], dictionary[kMaxPads], nCont = 0, padnumber = 0;
295 0 : Int_t *iter = &hvals[0], *citer = &content[0], *diter = &dictionary[0]; // iterators for preselection
296 : const Int_t threshold = 2;
297 0 : hvals[nPads] = -1; // termination for iterator
298 0 : do{
299 0 : if(*iter > threshold){
300 0 : *(citer++) = *iter;
301 0 : *(diter++) = padnumber;
302 0 : nCont++;
303 0 : }
304 0 : padnumber++;
305 0 : }while(*(++iter) != -1);
306 0 : TMath::Sort(nCont, content, indices);
307 :
308 : Int_t col, row, lower, lower1, upper, upper1;
309 0 : for(Int_t ib = 0; ib < nCont; ib++){
310 0 : if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
311 0 : AliDebug(1, Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack));
312 0 : break;
313 : }
314 : // Positions
315 0 : row = dictionary[indices[ib]]/nCols;
316 0 : col = dictionary[indices[ib]]%nCols;
317 : // here will be the threshold condition:
318 0 : if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
319 : // if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
320 : // break; // number of clusters below threshold: break;
321 : // }
322 : // passing: Mark the neighbors
323 0 : lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
324 0 : lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
325 0 : for(Int_t ic = lower; ic < upper; ++ic)
326 0 : for(Int_t ir = lower1; ir < upper1; ++ir){
327 0 : if(ic == col && ir == row) continue;
328 0 : mask[ic] |= (1 << ir);
329 0 : }
330 : // Storing the position in an array
331 : // testing for neigboring
332 : cogv = 0;
333 : norm = 0;
334 0 : lower = TMath::Max(col - 1, 0);
335 0 : upper = TMath::Min(col + 2, nCols);
336 0 : for(Int_t inb = lower; inb < upper; ++inb){
337 0 : cogv += yLengths[inb] * histogram[row][inb];
338 0 : norm += histogram[row][inb];
339 : }
340 0 : cogy[row][col] = cogv / norm;
341 : cogv = 0; norm = 0;
342 0 : lower = TMath::Max(row - 1, 0);
343 0 : upper = TMath::Min(row + 2, nRows);
344 0 : for(Int_t inb = lower; inb < upper; ++inb){
345 0 : cogv += zLengths[inb] * histogram[inb][col];
346 0 : norm += histogram[inb][col];
347 : }
348 0 : cogz[row][col] = Float_t(cogv) / norm;
349 : // passed the filter
350 0 : cand[nCandidates] = row*nCols + col; // store the position of a passig candidate into an Array
351 0 : sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
352 : // Analysis output
353 0 : nCandidates++;
354 0 : }
355 0 : if(!nCandidates) return kFALSE;
356 :
357 0 : Float_t pos[3], sig[2];
358 0 : Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
359 :
360 0 : new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
361 0 : fakeLayer->SetReconstructor(rec);
362 0 : fakeLayer->SetNRows(nRows);
363 0 : fakeLayer->SetOwner(kFALSE);
364 0 : if(nCandidates){
365 : UInt_t fakeIndex = 0;
366 0 : for(Int_t ican = 0; ican < nCandidates; ican++){
367 0 : row = cand[ican] / nCols;
368 0 : col = cand[ican] % nCols;
369 : //temporary
370 : Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
371 0 : for(int itb=0; itb<AliTRDseedV1::kNtb; itb++){
372 0 : if(!(nClusters = fTB[itb].GetNClusters())) continue;
373 0 : for(Int_t incl = 0; incl < nClusters; incl++){
374 0 : c = fTB[itb].GetCluster(incl);
375 0 : if(c->GetPadRow() != row) continue;
376 0 : if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
377 0 : x += c->GetX();
378 0 : y += c->GetY();
379 0 : z += c->GetZ();
380 0 : n++;
381 0 : }
382 0 : }
383 0 : if(!n) continue;
384 0 : pos[0] = x/n;
385 0 : pos[1] = y/n;
386 0 : pos[2] = z/n;
387 0 : sig[0] = .02;
388 0 : sig[1] = sigcands[ican];
389 0 : fakeLayer->InsertCluster(new AliTRDcluster(fDetector, 0., pos, sig, NULL, 3, signal, col, row, 0, 0, 0., 0), fakeIndex++);
390 0 : }
391 0 : }
392 0 : fakeLayer->BuildIndices();
393 : //fakeLayer->Print();
394 :
395 0 : if(rec->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kTracker) >= 3){
396 : //TMatrixD hist(nRows, nCols);
397 : //for(Int_t i = 0; i < nRows; i++)
398 : // for(Int_t j = 0; j < nCols; j++)
399 : // hist(i,j) = histogram[i][j];
400 0 : TTreeSRedirector &cstreamer = *rec->GetDebugStream(AliTRDrecoParam::kTracker);
401 0 : cstreamer << "GetSeedingLayer"
402 0 : << "layer=" << layer
403 0 : << "ymin=" << ymin
404 0 : << "ymax=" << ymax
405 0 : << "zmin=" << zmin
406 0 : << "zmax=" << zmax
407 0 : << "L.=" << fakeLayer
408 : //<< "Histogram.=" << &hist
409 0 : << "\n";
410 0 : }
411 :
412 : return kTRUE;
413 0 : }
414 :
415 :
416 : //_______________________________________________________
417 : void AliTRDtrackingChamber::Print(Option_t *opt) const
418 : {
419 : // Print the chamber status
420 0 : if(!GetNClusters()) return;
421 0 : AliInfo(Form("fDetector = %d", fDetector));
422 0 : AliInfo(Form("fX0 = %7.3f", fX0));
423 0 : const AliTRDchamberTimeBin *itb = &fTB[0];
424 0 : for(Int_t jtb=0; jtb<AliTRDseedV1::kNtb; jtb++, itb++) (*itb).Print(opt);
425 0 : }
426 :
427 :
428 : //_______________________________________________________
429 : void AliTRDtrackingChamber::Update()
430 : {
431 : // Steer purging of used and shared clusters
432 :
433 0 : AliTRDchamberTimeBin *jtb = &fTB[0];
434 0 : for(Int_t itb=AliTRDseedV1::kNtb; itb--; ++jtb){
435 0 : if(!(Int_t(*jtb))) continue;
436 0 : (*jtb).BuildIndices();
437 0 : }
438 0 : }
439 :
|
