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