Add some more parameters for stand alone cosmics reconstruction
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackingChamber.cxx
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
16 /* $Id: AliTRDtrackingChamber.cxx 23810 2008-02-08 09:00:27Z hristov $ */\r
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
41 ClassImp(AliTRDtrackingChamber)\r
42 \r
43 //_______________________________________________________\r
44 AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :\r
45         fDetector(det)\r
46         ,fX0(0.)\r
47 {}  \r
48 \r
49 //_______________________________________________________\r
50 void 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
56 void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)\r
57 {\r
58         fTB[c->GetLocalTimeBin()].InsertCluster(c, index);\r
59 }\r
60 \r
61 //_______________________________________________________\r
62 Bool_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
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
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
95 Int_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
107 Double_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
159 Bool_t AliTRDtrackingChamber::GetSeedingLayer(AliTRDchamberTimeBin *&fakeLayer, AliTRDgeometry *geo)\r
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
169         Int_t timeBinMin = AliTRDReconstructor::RecoParam()->GetNumberOfPresamples();\r
170         Int_t timeBinMax = AliTRDReconstructor::RecoParam()->GetNumberOfPostsamples();\r
171 \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
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
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
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
246         Int_t threshold = Int_t(maximum * AliTRDReconstructor::RecoParam()->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
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
294         if(!nCandidates) return kFALSE;\r
295         \r
296         Float_t pos[3], sig[2];\r
297         Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));\r
298         \r
299   new(fakeLayer) AliTRDchamberTimeBin(layer, stack, sector, z0, zl);\r
300         AliTRDcluster *cluster = 0x0;\r
301         if(nCandidates){\r
302                 UInt_t fakeIndex = 0;\r
303                 for(Int_t ican = 0; ican < nCandidates; ican++){\r
304                         row = cand[ican] / nCols;\r
305                         col = cand[ican] % nCols;\r
306                         //temporary\r
307                         Int_t n = 0; Double_t x = 0., y = 0., z = 0.;\r
308                         for(int itb=0; itb<kNTimeBins; itb++){\r
309                                 if(!(nClusters = fTB[itb].GetNClusters())) continue;\r
310                                 for(Int_t incl = 0; incl < nClusters; incl++){\r
311                                         c = fTB[itb].GetCluster(incl);  \r
312                                         if(c->GetPadRow() != row) continue;\r
313                                         if(TMath::Abs(c->GetPadCol() - col) > 2) continue;\r
314                                         x += c->GetX();\r
315                                         y += c->GetY();\r
316                                         z += c->GetZ();\r
317                                         n++;\r
318                                 }\r
319                         }\r
320                         pos[0] = x/n;\r
321                         pos[1] = y/n;\r
322                         pos[2] = z/n;\r
323                         sig[0] = .02;\r
324                         sig[1] = sigcands[ican];\r
325                         cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);\r
326                         fakeLayer->InsertCluster(cluster, fakeIndex++);\r
327                 }\r
328         }\r
329         fakeLayer->SetNRows(nRows);\r
330         fakeLayer->SetOwner();\r
331         fakeLayer->BuildIndices();\r
332         //fakeLayer->PrintClusters();\r
333         \r
334         if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 3){\r
335                 //TMatrixD hist(nRows, nCols);\r
336                 //for(Int_t i = 0; i < nRows; i++)\r
337                 //      for(Int_t j = 0; j < nCols; j++)\r
338                 //              hist(i,j) = histogram[i][j];\r
339                 TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();\r
340                 cstreamer << "GetSeedingLayer"\r
341                 << "layer="      << layer\r
342                 << "ymin="       << ymin\r
343                 << "ymax="       << ymax\r
344                 << "zmin="       << zmin\r
345                 << "zmax="       << zmax\r
346                 << "L.="         << fakeLayer\r
347                 //<< "Histogram.=" << &hist\r
348                 << "\n";\r
349         }\r
350         \r
351         return kTRUE;\r
352 }\r
353 \r