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