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