]> 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
Update database entry
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackingChamber.cxx
CommitLineData
d20df6fc 1/**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\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
15\r
4e459a9d 16/* $Id: AliTRDtrackingChamber.cxx 23810 2008-02-08 09:00:27Z hristov $ */\r
d20df6fc 17\r
18////////////////////////////////////////////////////////////////////////////\r
19// //\r
20// Tracking in one chamber //\r
21// //\r
22// Authors: //\r
23// Alex Bercuci <A.Bercuci@gsi.de> //\r
24// Markus Fasel <M.Fasel@gsi.de> //\r
25// //\r
26////////////////////////////////////////////////////////////////////////////\r
27\r
28#include "AliTRDtrackingChamber.h"\r
29\r
30#include "TMath.h"\r
31#include "TMatrixTBase.h"\r
32#include <TTreeStream.h>\r
33\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
40\r
41ClassImp(AliTRDtrackingChamber)\r
42\r
43//_______________________________________________________\r
44AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :\r
45 fDetector(det)\r
46 ,fX0(0.)\r
47{} \r
48\r
49//_______________________________________________________\r
50void AliTRDtrackingChamber::Clear(const Option_t *opt)\r
51{\r
52 for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);\r
53}\r
54\r
55//_______________________________________________________\r
56void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)\r
57{\r
3a039a31 58 fTB[c->GetPadTime()].InsertCluster(c, index);\r
d20df6fc 59}\r
60\r
61//_______________________________________________________\r
62Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo)\r
63{\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
67//\r
053767a4 68 Int_t stack = geo->GetStack(fDetector);\r
69 Int_t layer = geo->GetLayer(fDetector);\r
70 AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);\r
d20df6fc 71 Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();\r
053767a4 72 Double_t z0 = geo->GetRow0(layer, stack, 0) - zl;\r
d20df6fc 73 Int_t nrows = pp->GetNrows();\r
74 \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
81 index[jtb++] = itb;\r
82 } \r
83 if(jtb<2) return kFALSE;\r
84 \r
85 \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
91 return kTRUE;\r
92}\r
93 \r
94//_______________________________________________________ \r
95Int_t AliTRDtrackingChamber::GetNClusters() const\r
96{\r
97// Returns number of clusters in chamber\r
98//\r
99 Int_t n = 0;\r
100 for(Int_t itb=0; itb<kNTimeBins; itb++){ \r
101 n += Int_t(fTB[itb]);\r
102 }\r
103 return n; \r
104} \r
105\r
106//_______________________________________________________\r
107Double_t AliTRDtrackingChamber::GetQuality()\r
108{\r
109 //\r
110 // Calculate chamber quality for seeding.\r
111 // \r
112 //\r
113 // Parameters :\r
114 // layers : Array of propagation layers for this plane.\r
115 //\r
116 // Output :\r
117 // plane quality factor for seeding\r
118 // \r
119 // Detailed description\r
120 //\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
126 //\r
127\r
128 Int_t ncl = 0;\r
129 Int_t nused = 0;\r
130 Int_t nClLayer;\r
131 for(int itb=0; itb<kNTimeBins; itb++){\r
132 if(!(nClLayer = fTB[itb].GetNClusters())) continue;\r
133 ncl += nClLayer;\r
134 for(Int_t incl = 0; incl < nClLayer; incl++){\r
135 if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;\r
136 }\r
137 }\r
138 \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
146\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
151// \r
152// printf("ncl_dev(%f) a(%f)\n", ncl_dev, a);\r
153// return quality*TMath::Exp(-a);\r
154\r
155}\r
156\r
157\r
158//_______________________________________________________\r
3a039a31 159Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry *geo, const AliTRDReconstructor *rec)\r
d20df6fc 160{\r
161 //\r
162 // Creates a seeding layer\r
163 //\r
164 \r
165 // constants\r
166 const Int_t kMaxRows = 16;\r
167 const Int_t kMaxCols = 144;\r
168 const Int_t kMaxPads = 2304;\r
3a039a31 169 Int_t timeBinMin = rec->GetRecoParam()->GetNumberOfPresamples();\r
170 Int_t timeBinMax = rec->GetRecoParam()->GetNumberOfPostsamples();\r
4e459a9d 171\r
d20df6fc 172 // Get the geometrical data of the chamber\r
053767a4 173 Int_t layer = geo->GetLayer(fDetector);\r
174 Int_t stack = geo->GetStack(fDetector);\r
d20df6fc 175 Int_t sector= geo->GetSector(fDetector);\r
053767a4 176 AliTRDpadPlane *pp = geo->GetPadPlane(layer, stack);\r
d20df6fc 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
186 \r
187 // Fill the histogram\r
188 Int_t nClusters; \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
197 }\r
4e459a9d 198 for(Int_t iTime = timeBinMin; iTime < kNTimeBins-timeBinMax; iTime++){\r
199 if(!(nClusters = fTB[iTime].GetNClusters())) continue;\r
d20df6fc 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
206 }\r
207 }\r
208 \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
216 \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
222 }\r
223 for(Int_t icnt = 0; icnt < nRows; icnt++){\r
224 zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;\r
225 }\r
226\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
232 }\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
236 \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
244 // Set Threshold\r
245 Int_t maximum = hvals[indices[0]]; // best\r
3a039a31 246 Int_t threshold = Int_t(maximum * rec->GetRecoParam()->GetFindableClusters());\r
d20df6fc 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
251 break;\r
252 }\r
253 // Positions\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
260 } \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
268 }\r
269 // Storing the position in an array\r
270 // testing for neigboring\r
271 cogv = 0;\r
272 norm = 0;\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
278 }\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
286 }\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
291 // Analysis output\r
292 nCandidates++;\r
293 }\r
d611c74f 294 if(!nCandidates) return kFALSE;\r
d20df6fc 295 \r
296 Float_t pos[3], sig[2];\r
297 Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));\r
d611c74f 298 \r
299 new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);\r
43d6ad34 300 fakeLayer->SetReconstructor(rec);\r
d20df6fc 301 AliTRDcluster *cluster = 0x0;\r
302 if(nCandidates){\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
307 //temporary\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
315 x += c->GetX();\r
316 y += c->GetY();\r
317 z += c->GetZ();\r
318 n++;\r
319 }\r
320 }\r
321 pos[0] = x/n;\r
322 pos[1] = y/n;\r
323 pos[2] = z/n;\r
324 sig[0] = .02;\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
328 }\r
329 }\r
330 fakeLayer->SetNRows(nRows);\r
331 fakeLayer->SetOwner();\r
332 fakeLayer->BuildIndices();\r
333 //fakeLayer->PrintClusters();\r
334 \r
3a039a31 335 if(rec->GetStreamLevel(AliTRDReconstructor::kTracker) >= 3){\r
d20df6fc 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
053767a4 342 << "layer=" << layer\r
d20df6fc 343 << "ymin=" << ymin\r
344 << "ymax=" << ymax\r
345 << "zmin=" << zmin\r
346 << "zmax=" << zmax\r
347 << "L.=" << fakeLayer\r
348 //<< "Histogram.=" << &hist\r
349 << "\n";\r
350 }\r
351 \r
d611c74f 352 return kTRUE;\r
d20df6fc 353}\r
354\r