1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
16 /* $Id: AliTRDtrackingChamber.cxx 23810 2008-02-08 09:00:27Z hristov $ */
\r
18 ////////////////////////////////////////////////////////////////////////////
\r
20 // Tracking in one chamber //
\r
23 // Alex Bercuci <A.Bercuci@gsi.de> //
\r
24 // Markus Fasel <M.Fasel@gsi.de> //
\r
26 ////////////////////////////////////////////////////////////////////////////
\r
28 #include "AliTRDtrackingChamber.h"
\r
31 #include "TMatrixTBase.h"
\r
32 #include <TTreeStream.h>
\r
34 #include "AliTRDReconstructor.h"
\r
35 #include "AliTRDrecoParam.h"
\r
36 #include "AliTRDtrackerV1.h"
\r
37 #include "AliTRDgeometry.h"
\r
38 #include "AliTRDpadPlane.h"
\r
39 #include "AliTRDcalibDB.h"
\r
41 ClassImp(AliTRDtrackingChamber)
\r
43 //_______________________________________________________
\r
44 AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :
\r
49 //_______________________________________________________
\r
50 void AliTRDtrackingChamber::Clear(const Option_t *opt)
\r
52 for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);
\r
55 //_______________________________________________________
\r
56 void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)
\r
58 fTB[c->GetPadTime()].InsertCluster(c, index);
\r
61 //_______________________________________________________
\r
62 Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo)
\r
64 // Init chamber and all time bins (AliTRDchamberTimeBin)
\r
65 // Calculates radial position of the chamber based on
\r
66 // radial positions of the time bins (calibration/alignment aware)
\r
68 Int_t stack = geo->GetStack(fDetector);
\r
69 Int_t layer = geo->GetLayer(fDetector);
\r
70 AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
\r
71 Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
\r
72 Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;
\r
73 Int_t nrows = pp->GetNrows();
\r
75 Int_t index[50], jtb = 0;
\r
76 for(Int_t itb=0; itb<kNTimeBins; itb++){
\r
77 if(!fTB[itb]) continue;
\r
78 fTB[itb].SetRange(z0, zl);
\r
79 fTB[itb].SetNRows(nrows);
\r
80 fTB[itb].BuildIndices();
\r
83 if(jtb<2) return kFALSE;
\r
86 // ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
\r
87 Double_t x0 = fTB[index[0]].GetX();
\r
88 Double_t x1 = fTB[index[1]].GetX();
\r
89 Double_t dx = (x0 - x1)/(index[1] - index[0]);
\r
90 fX0 = x0 + dx*(index[0] - (Int_t)AliTRDcalibDB::Instance()->GetT0Average(fDetector));
\r
94 //_______________________________________________________
\r
95 Int_t AliTRDtrackingChamber::GetNClusters() const
\r
97 // Returns number of clusters in chamber
\r
100 for(Int_t itb=0; itb<kNTimeBins; itb++){
\r
101 n += Int_t(fTB[itb]);
\r
106 //_______________________________________________________
\r
107 Double_t AliTRDtrackingChamber::GetQuality()
\r
110 // Calculate chamber quality for seeding.
\r
114 // layers : Array of propagation layers for this plane.
\r
117 // plane quality factor for seeding
\r
119 // Detailed description
\r
121 // The quality of the plane for seeding is higher if:
\r
122 // 1. the average timebin population is closer to an integer number
\r
123 // 2. the distribution of clusters/timebin is closer to a uniform distribution.
\r
124 // - the slope of the first derivative of a parabolic fit is small or
\r
125 // - the slope of a linear fit is small
\r
131 for(int itb=0; itb<kNTimeBins; itb++){
\r
132 if(!(nClLayer = fTB[itb].GetNClusters())) continue;
\r
134 for(Int_t incl = 0; incl < nClLayer; incl++){
\r
135 if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
\r
139 // calculate the deviation of the mean number of clusters from the
\r
140 // closest integer values
\r
141 Float_t nclMed = float(ncl-nused)/AliTRDtrackerV1::GetNTimeBins();
\r
142 Int_t ncli = Int_t(nclMed);
\r
143 Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
\r
144 nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
\r
145 return TMath::Exp(-5.*TMath::Abs(nclDev));
\r
147 // // get slope of the derivative
\r
148 // if(!fitter.Eval()) return quality;
\r
149 // fitter.PrintResults(3);
\r
150 // Double_t a = fitter.GetParameter(1);
\r
152 // printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);
\r
153 // return quality*TMath::Exp(-a);
\r
158 //_______________________________________________________
\r
159 Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry *geo, const AliTRDReconstructor *rec)
\r
162 // Creates a seeding layer
\r
166 const Int_t kMaxRows = 16;
\r
167 const Int_t kMaxCols = 144;
\r
168 const Int_t kMaxPads = 2304;
\r
169 Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();
\r
170 Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();
\r
172 // Get the geometrical data of the chamber
\r
173 Int_t layer = geo->GetLayer(fDetector);
\r
174 Int_t stack = geo->GetStack(fDetector);
\r
175 Int_t sector= geo->GetSector(fDetector);
\r
176 AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);
\r
177 Int_t nCols = pp->GetNcols();
\r
178 Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
\r
179 Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
\r
180 Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
\r
181 Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
\r
182 Float_t z0 = -1., zl = -1.;
\r
183 Int_t nRows = pp->GetNrows();
\r
184 Float_t binlength = (ymax - ymin)/nCols;
\r
185 //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
\r
187 // Fill the histogram
\r
189 Int_t *histogram[kMaxRows]; // 2D-Histogram
\r
190 Int_t hvals[kMaxPads]; memset(hvals, 0, sizeof(Int_t)*kMaxPads);
\r
191 Float_t *sigmas[kMaxRows];
\r
192 Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads);
\r
193 AliTRDcluster *c = 0x0;
\r
194 for(Int_t irs = 0; irs < kMaxRows; irs++){
\r
195 histogram[irs] = &hvals[irs*kMaxCols];
\r
196 sigmas[irs] = &svals[irs*kMaxCols];
\r
198 for(Int_t iTime = timeBinMin; iTime < kNTimeBins-timeBinMax; iTime++){
\r
199 if(!(nClusters = fTB[iTime].GetNClusters())) continue;
\r
200 z0 = fTB[iTime].GetZ0();
\r
201 zl = fTB[iTime].GetDZ0();
\r
202 for(Int_t incl = 0; incl < nClusters; incl++){
\r
203 c = fTB[iTime].GetCluster(incl);
\r
204 histogram[c->GetPadRow()][c->GetPadCol()]++;
\r
205 sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
\r
209 // Now I have everything in the histogram, do the selection
\r
210 //Int_t nPads = nCols * nRows;
\r
211 // This is what we are interested in: The center of gravity of the best candidates
\r
212 Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
\r
213 Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
\r
214 Float_t *cogy[kMaxRows];
\r
215 Float_t *cogz[kMaxRows];
\r
217 // Lookup-Table storing coordinates according to the bins
\r
218 Float_t yLengths[kMaxCols];
\r
219 Float_t zLengths[kMaxRows];
\r
220 for(Int_t icnt = 0; icnt < nCols; icnt++){
\r
221 yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
\r
223 for(Int_t icnt = 0; icnt < nRows; icnt++){
\r
224 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
\r
227 // A bitfield is used to mask the pads as usable
\r
228 Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
\r
229 for(UChar_t icount = 0; icount < nRows; icount++){
\r
230 cogy[icount] = &cogyvals[icount*kMaxCols];
\r
231 cogz[icount] = &cogzvals[icount*kMaxCols];
\r
233 // In this array the array position of the best candidates will be stored
\r
234 Int_t cand[AliTRDtrackerV1::kMaxTracksStack];
\r
235 Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
\r
237 // helper variables
\r
238 Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
\r
239 Int_t nCandidates = 0;
\r
240 Float_t norm, cogv;
\r
241 // histogram filled -> Select best bins
\r
242 Int_t nPads = nCols * nRows;
\r
243 TMath::Sort(nPads, hvals, indices); // bins storing a 0 should not matter
\r
245 Int_t maximum = hvals[indices[0]]; // best
\r
246 Int_t threshold = Int_t(maximum * rec->GetRecoParam()->GetFindableClusters());
\r
247 Int_t col, row, lower, lower1, upper, upper1;
\r
248 for(Int_t ib = 0; ib < nPads; ib++){
\r
249 if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
\r
250 printf("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::kMaxTracksStack);
\r
254 row = indices[ib]/nCols;
\r
255 col = indices[ib]%nCols;
\r
256 // here will be the threshold condition:
\r
257 if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue
\r
258 if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed
\r
259 break; // number of clusters below threshold: break;
\r
261 // passing: Mark the neighbors
\r
262 lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols);
\r
263 lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
\r
264 for(Int_t ic = lower; ic < upper; ++ic)
\r
265 for(Int_t ir = lower1; ir < upper1; ++ir){
\r
266 if(ic == col && ir == row) continue;
\r
267 mask[ic] |= (1 << ir);
\r
269 // Storing the position in an array
\r
270 // testing for neigboring
\r
273 lower = TMath::Max(col - 1, 0);
\r
274 upper = TMath::Min(col + 2, nCols);
\r
275 for(Int_t inb = lower; inb < upper; ++inb){
\r
276 cogv += yLengths[inb] * histogram[row][inb];
\r
277 norm += histogram[row][inb];
\r
279 cogy[row][col] = cogv / norm;
\r
280 cogv = 0; norm = 0;
\r
281 lower = TMath::Max(row - 1, 0);
\r
282 upper = TMath::Min(row + 2, nRows);
\r
283 for(Int_t inb = lower; inb < upper; ++inb){
\r
284 cogv += zLengths[inb] * histogram[inb][col];
\r
285 norm += histogram[inb][col];
\r
287 cogz[row][col] = Float_t(cogv) / norm;
\r
288 // passed the filter
\r
289 cand[nCandidates] = row*nCols + col; // store the position of a passig candidate into an Array
\r
290 sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
\r
294 if(!nCandidates) return kFALSE;
\r
296 Float_t pos[3], sig[2];
\r
297 Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
\r
299 new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);
\r
300 fakeLayer->SetReconstructor(rec);
\r
301 AliTRDcluster *cluster = 0x0;
\r
303 UInt_t fakeIndex = 0;
\r
304 for(Int_t ican = 0; ican < nCandidates; ican++){
\r
305 row = cand[ican] / nCols;
\r
306 col = cand[ican] % nCols;
\r
308 Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
\r
309 for(int itb=0; itb<kNTimeBins; itb++){
\r
310 if(!(nClusters = fTB[itb].GetNClusters())) continue;
\r
311 for(Int_t incl = 0; incl < nClusters; incl++){
\r
312 c = fTB[itb].GetCluster(incl);
\r
313 if(c->GetPadRow() != row) continue;
\r
314 if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
\r
325 sig[1] = sigcands[ican];
\r
326 cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);
\r
327 fakeLayer->InsertCluster(cluster, fakeIndex++);
\r
330 fakeLayer->SetNRows(nRows);
\r
331 fakeLayer->SetOwner();
\r
332 fakeLayer->BuildIndices();
\r
333 //fakeLayer->PrintClusters();
\r
335 if(rec->GetStreamLevel(AliTRDReconstructor::kTracker) >= 3){
\r
336 //TMatrixD hist(nRows, nCols);
\r
337 //for(Int_t i = 0; i < nRows; i++)
\r
338 // for(Int_t j = 0; j < nCols; j++)
\r
339 // hist(i,j) = histogram[i][j];
\r
340 TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
\r
341 cstreamer << "GetSeedingLayer"
\r
342 << "layer=" << layer
\r
347 << "L.=" << fakeLayer
\r
348 //<< "Histogram.=" << &hist
\r