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