Fix bug in the definition of number of clusters/tracklet
[u/mrichter/AliRoot.git] / TRD / AliTRDseedV1.h
1 #ifndef ALITRDSEEDV1_H
2 #define ALITRDSEEDV1_H
3 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * See cxx source for full Copyright notice                               */
5
6 /* $Id$ */
7
8 ////////////////////////////////////////////////////////////////////////////
9 //                                                                        //
10 // \class AliTRDseedV1
11 // \brief The TRD offline tracklet
12 // \author Alexandru Bercuci
13 //                                                                        //
14 ////////////////////////////////////////////////////////////////////////////
15
16 #ifndef ALITRDTRACKLETBASE_H
17 #include "AliTRDtrackletBase.h"
18 #endif
19
20 #ifndef ROOT_TMath
21 #include "TMath.h"
22 #endif
23
24 #ifndef ALITRDGEOMETRY_H
25 #include "AliTRDgeometry.h"
26 #endif
27
28 #ifndef ALIPID_H
29 #include "AliPID.h"
30 #endif
31
32 #ifndef ALIRIEMAN_H
33 #include "AliRieman.h"
34 #endif
35
36 #ifndef ALITRDCLUSTER_H 
37 #include "AliTRDcluster.h"
38 #endif
39
40 #include "AliTRDReconstructor.h"
41
42 class TTreeSRedirector;
43 class TLinearFitter;
44
45 class AliRieman;
46
47 class AliTRDtrackingChamber;
48 class AliTRDtrackV1;
49 class AliTRDpadPlane;
50 class AliTRDseedV1 : public AliTRDtrackletBase
51 {
52   friend class AliHLTTRDTracklet;
53
54 public:
55   enum ETRDtrackletBuffers {    
56     kNtb       = 31     // max clusters/pad row
57    ,kNclusters = 2*kNtb // max number of clusters/tracklet
58    ,kNslices   = 10     // max dEdx slices
59   };
60
61   // bits from 0-13 are reserved by ROOT (see TObject.h)
62   enum ETRDtrackletStatus {
63     kOwner      = BIT(14) // owner of its clusters
64    ,kRowCross   = BIT(15) // pad row cross tracklet
65    ,kPID        = BIT(16) // PID contributor
66    ,kCalib      = BIT(17) // calibrated tracklet
67    ,kKink       = BIT(18) // kink prolongation tracklet
68    ,kStandAlone = BIT(19) // tracklet build during stand alone track finding
69   };
70   enum ETRDtrackletError {
71     kAttachClFound = 1  // not enough clusters found
72    ,kAttachRow          // found row < 0
73    ,kAttachMultipleCl   // multiple clusters attached to time bin
74    ,kAttachClAttach     // not enough clusters attached
75   };
76
77   AliTRDseedV1(Int_t det = -1);
78   ~AliTRDseedV1();
79   AliTRDseedV1(const AliTRDseedV1 &ref);
80   AliTRDseedV1& operator=(const AliTRDseedV1 &ref);
81
82   Bool_t    AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt = kFALSE);
83   void      Bootstrap(const AliTRDReconstructor *rec);
84   void      Calibrate();
85   void      CookdEdx(Int_t nslices);
86   void      CookLabels();
87   Bool_t    CookPID();
88   Bool_t    Fit(Bool_t tilt=kFALSE, Bool_t zcorr=kFALSE);
89   Bool_t    Init(AliTRDtrackV1 *track);
90   inline void      Init(const AliRieman *fit);
91   Bool_t    IsEqual(const TObject *inTracklet) const;
92   Bool_t    IsCalibrated() const     { return TestBit(kCalib);}
93   Bool_t    IsOwner() const          { return TestBit(kOwner);}
94   Bool_t    IsKink() const           { return TestBit(kKink);}
95   Bool_t    HasPID() const           { return TestBit(kPID);}
96   Bool_t    IsOK() const             { return GetN() > 4 && GetNUsed() < 4;}
97   Bool_t    IsRowCross() const       { return TestBit(kRowCross);}
98   Bool_t    IsUsable(Int_t i) const  { return fClusters[i] && !fClusters[i]->IsUsed();}
99   Bool_t    IsStandAlone() const     { return TestBit(kStandAlone);}
100
101   Float_t   GetC() const             { return fC; }
102   Float_t   GetChi2() const          { return fChi2; }
103   inline Float_t   GetChi2Z() const;
104   inline Float_t   GetChi2Y() const;
105   inline Float_t   GetChi2Phi() const;
106   void      GetCovAt(Double_t x, Double_t *cov) const;
107   void      GetCovXY(Double_t *cov) const { memcpy(cov, &fCov[0], 3*sizeof(Double_t));}
108   void      GetCovRef(Double_t *cov) const { memcpy(cov, &fRefCov, 7*sizeof(Double_t));}
109   static Double_t GetCovSqrt(const Double_t * const c, Double_t *d);
110   static Double_t GetCovInv(const Double_t * const c, Double_t *d);
111   UChar_t   GetErrorMsg() const      { return fErrorMsg;}
112   Float_t   GetdX() const            { return fdX;}
113   const Float_t*  GetdEdx() const    { return &fdEdx[0];}
114   Float_t   GetdQdl(Int_t ic, Float_t *dx=NULL) const;
115   Float_t   GetdYdX() const          { return fYfit[1]; } 
116   Float_t   GetdZdX() const          { return fZref[1]; }
117   Int_t     GetdY() const            { return Int_t(GetY()/0.014);}
118   Int_t     GetDetector() const      { return fDet;}
119   void      GetCalibParam(Float_t &exb, Float_t &vd, Float_t &t0, Float_t &s2, Float_t &dl, Float_t &dt) const    { 
120               exb = fExB; vd = fVD; t0 = fT0; s2 = fS2PRF; dl = fDiffL; dt = fDiffT;}
121   AliTRDcluster*  GetClusters(Int_t i) const               { return i<0 || i>=kNclusters ? NULL: fClusters[i];}
122   static TLinearFitter*  GetFitterY();
123   static TLinearFitter*  GetFitterZ();
124   Int_t     GetIndexes(Int_t i) const{ return i<0 || i>=kNclusters ? -1 : fIndexes[i];}
125   Int_t     GetLabels(Int_t i) const { return fLabels[i];}  
126   Float_t   GetMomentum(Float_t *err = NULL) const;
127   Int_t     GetN() const             { return (Int_t)fN&0x1f;}
128   Int_t     GetN2() const            { return GetN();}
129   Int_t     GetNUsed() const         { return Int_t((fN>>5)&0x1f);}
130   Int_t     GetNShared() const       { return Int_t((fN>>10)&0x1f);}
131   Float_t   GetQuality(Bool_t kZcorr) const;
132   Float_t   GetPadLength() const     { return fPad[0];}
133   Float_t   GetPadWidth() const      { return fPad[1];}
134   Int_t     GetPlane() const         { return AliTRDgeometry::GetLayer(fDet);    }
135
136   Float_t*  GetProbability(Bool_t force=kFALSE);
137   Float_t   GetPt() const            { return fPt; }
138   inline Double_t  GetPID(Int_t is=-1) const;
139   Float_t   GetS2Y() const           { return fS2Y;}
140   Float_t   GetS2Z() const           { return fS2Z;}
141   Float_t   GetSigmaY() const        { return fS2Y > 0. ? TMath::Sqrt(fS2Y) : 0.2;}
142   Float_t   GetSnp() const           { return fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);}
143   Float_t   GetTgl() const           { return fZref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);}
144   Float_t   GetTilt() const          { return fPad[2];}
145   UInt_t    GetTrackletWord() const  { return 0;}
146   UShort_t  GetVolumeId() const;
147   Float_t   GetX0() const            { return fX0;}
148   Float_t   GetX() const             { return fX0 - fX;}
149   Float_t   GetY() const             { return fYfit[0] - fYfit[1] * fX;}
150   Double_t  GetYat(Double_t x) const { return fYfit[0] - fYfit[1] * (fX0-x);}
151   Float_t   GetYfit(Int_t id) const  { return fYfit[id];}
152   Float_t   GetYref(Int_t id) const  { return fYref[id];}
153   Float_t   GetZ() const             { return fZfit[0] - fZfit[1] * fX;}
154   Double_t  GetZat(Double_t x) const { return fZfit[0] - fZfit[1] * (fX0-x);}
155   Float_t   GetZfit(Int_t id) const  { return fZfit[id];}
156   Float_t   GetZref(Int_t id) const  { return fZref[id];}
157   Int_t     GetYbin() const          { return Int_t(GetY()/0.016);}
158   Int_t     GetZbin() const          { return Int_t(GetZ()/fPad[0]);}
159
160   inline AliTRDcluster* NextCluster();
161   inline AliTRDcluster* PrevCluster();
162   void      Print(Option_t *o = "") const;
163   inline void ResetClusterIter(Bool_t forward = kTRUE);
164   void      Reset();
165
166   void      SetC(Float_t c)          { fC = c;}
167   void      SetChi2(Float_t chi2)    { fChi2 = chi2;}
168   inline void SetCovRef(const Double_t *cov);
169   void      SetErrorMsg(Int_t err)   { fErrorMsg = err;}
170   void      SetIndexes(Int_t i, Int_t idx) { fIndexes[i]  = idx; }
171   void      SetLabels(Int_t *lbls)   { memcpy(fLabels, lbls, 3*sizeof(Int_t)); }
172   void      SetKink(Bool_t k = kTRUE){ SetBit(kKink, k);}
173   void      SetPID(Bool_t k = kTRUE) { SetBit(kPID, k);}
174   void      SetStandAlone(Bool_t st) { SetBit(kStandAlone, st); }
175   void      SetPt(Double_t pt)       { fPt = pt;}
176   void      SetOwner();
177   void      SetPadPlane(AliTRDpadPlane *p);
178   void      SetPadLength(Float_t l)  { fPad[0] = l;}
179   void      SetPadWidth(Float_t w)   { fPad[1] = w;}
180   void      SetTilt(Float_t tilt)    { fPad[2] = tilt; }
181   void      SetDetector(Int_t d)     { fDet = d;  }
182   void      SetDX(Float_t inDX)      { fdX = inDX;}
183   void      SetReconstructor(const AliTRDReconstructor *rec) {fkReconstructor = rec;}
184   void      SetX0(Float_t x0)        { fX0 = x0; }
185   void      SetYref(Int_t i, Float_t y) { fYref[i]     = y;}
186   void      SetZref(Int_t i, Float_t z) { fZref[i]     = z;}
187 //   void      SetUsabilityMap(Long_t um)  { fUsable = um; }
188   void      Update(const AliTRDtrackV1* trk);
189   void      UpdateUsed();
190   void      UseClusters();
191
192 protected:
193   void        Copy(TObject &ref) const;
194
195 private:
196   inline void SetN(Int_t n);
197   inline void SetNUsed(Int_t n);
198   inline void SetNShared(Int_t n);
199
200   const AliTRDReconstructor *fkReconstructor;//! local reconstructor
201   AliTRDcluster  **fClusterIter;            //! clusters iterator
202   Int_t            fIndexes[kNclusters];    //! Indexes
203   Float_t          fExB;                    //! tg(a_L) @ tracklet location
204   Float_t          fVD;                     //! drift velocity @ tracklet location
205   Float_t          fT0;                     //! time 0 @ tracklet location
206   Float_t          fS2PRF;                  //! sigma^2 PRF for xd->0 and phi=a_L 
207   Float_t          fDiffL;                  //! longitudinal diffusion coefficient
208   Float_t          fDiffT;                  //! transversal diffusion coefficient
209   Char_t           fClusterIdx;             //! clusters iterator
210   UChar_t          fErrorMsg;               // processing error
211   UInt_t           fN;                      // number of clusters attached/used/shared
212   Short_t          fDet;                    // TRD detector
213   AliTRDcluster   *fClusters[kNclusters];   // Clusters
214   Float_t          fPad[3];                 // local pad definition : length/width/tilt 
215   Float_t          fYref[2];                //  Reference y, dydx
216   Float_t          fZref[2];                //  Reference z, dz/dx
217   Float_t          fYfit[2];                //  Fit y, dy/dx
218   Float_t          fZfit[2];                //  Fit z
219   Float_t          fPt;                     //  Pt estimate @ tracklet [GeV/c]
220   Float_t          fdX;                     // length of time bin
221   Float_t          fX0;                     // anode wire position
222   Float_t          fX;                      // radial position of the tracklet
223   Float_t          fY;                      // r-phi position of the tracklet
224   Float_t          fZ;                      // z position of the tracklet
225   Float_t          fS2Y;                    // estimated resolution in the r-phi direction 
226   Float_t          fS2Z;                    // estimated resolution in the z direction 
227   Float_t          fC;                      // Curvature
228   Float_t          fChi2;                   // Global chi2  
229   Float_t          fdEdx[kNslices];         // dE/dx measurements for tracklet
230   Float_t          fProb[AliPID::kSPECIES]; // PID probabilities
231   Int_t            fLabels[3];              // most frequent MC labels and total number of different labels
232   Double_t         fRefCov[7];              // covariance matrix of the track in the yz plane + the rest of the diagonal elements
233   Double_t         fCov[3];                 // covariance matrix of the tracklet in the xy plane
234   static TLinearFitter   *fgFitterY;        // Linear Fitter for tracklet fit in xy-plane
235   static TLinearFitter   *fgFitterZ;        // Linear Fitter for tracklet fit in xz-plane
236
237   ClassDef(AliTRDseedV1, 8)                 // The offline TRD tracklet 
238 };
239
240 //____________________________________________________________
241 inline Float_t AliTRDseedV1::GetChi2Z() const
242 {
243   Double_t dz = fZref[0]-fZfit[0]; dz*=dz;
244   Double_t cov[3]; GetCovAt(fX, cov);
245   Double_t s2 = fRefCov[2]+cov[2];
246   return s2 > 0. ? dz/s2 : 0.; 
247 }
248
249 //____________________________________________________________
250 inline Float_t AliTRDseedV1::GetChi2Y() const
251 {
252   Double_t dy = fYref[0]-fYfit[0]; dy*=dy;
253   Double_t cov[3]; GetCovAt(fX, cov);
254   Double_t s2 = fRefCov[0]+cov[0];
255   return s2 > 0. ? dy/s2 : 0.; 
256 }
257
258 //____________________________________________________________
259 inline Float_t AliTRDseedV1::GetChi2Phi() const
260 {
261   Double_t dphi = fYref[1]-fYfit[1]; dphi*=dphi;
262   Double_t cov[3]; GetCovAt(fX, cov);
263   Double_t s2 = fRefCov[2]+cov[2];
264   return s2 > 0. ? dphi/s2 : 0.; 
265 }
266
267
268
269 //____________________________________________________________
270 inline Double_t AliTRDseedV1::GetPID(Int_t is) const
271 {
272   if(is<0) return fProb[AliPID::kElectron];
273   if(is<AliPID::kSPECIES) return fProb[is];
274   return 0.;
275 }
276
277 //____________________________________________________________
278 inline void AliTRDseedV1::Init(const AliRieman *rieman)
279 {
280   fZref[0] = rieman->GetZat(fX0);
281   fZref[1] = rieman->GetDZat(fX0);
282   fYref[0] = rieman->GetYat(fX0);
283   fYref[1] = rieman->GetDYat(fX0);
284   if(fkReconstructor && fkReconstructor->IsHLT()){
285     fRefCov[0] = 1;
286     fRefCov[2] = 10;
287   }else{
288     fRefCov[0] = rieman->GetErrY(fX0);
289     fRefCov[2] = rieman->GetErrZ(fX0);
290   }
291   fC       = rieman->GetC(); 
292   fChi2    = rieman->GetChi2();
293 }
294
295 //____________________________________________________________
296 inline AliTRDcluster* AliTRDseedV1::NextCluster()
297 {
298 // Mimic the usage of STL iterators.
299 // Forward iterator
300
301   fClusterIdx++; fClusterIter++;
302   while(fClusterIdx < kNclusters){
303     if(!(*fClusterIter)){ 
304       fClusterIdx++; 
305       fClusterIter++;
306       continue;
307     }
308     return *fClusterIter;
309   }
310   return NULL;
311 }
312
313 //____________________________________________________________
314 inline AliTRDcluster* AliTRDseedV1::PrevCluster()
315 {
316 // Mimic the usage of STL iterators.
317 // Backward iterator
318
319   fClusterIdx--; fClusterIter--;
320   while(fClusterIdx >= 0){
321     if(!(*fClusterIter)){ 
322       fClusterIdx--; 
323       fClusterIter--;
324       continue;
325     }
326     return *fClusterIter;
327   }
328   return NULL;
329 }
330
331 //____________________________________________________________
332 inline void AliTRDseedV1::ResetClusterIter(Bool_t forward) 
333 {
334 // Mimic the usage of STL iterators.
335 // Facilitate the usage of NextCluster for forward like 
336 // iterator (kTRUE) and PrevCluster for backward like iterator (kFALSE)
337
338   if(forward){
339     fClusterIter = &fClusters[0]; fClusterIter--; 
340     fClusterIdx=-1;
341   } else {
342     fClusterIter = &fClusters[kNclusters-1]; fClusterIter++; 
343     fClusterIdx=kNclusters;
344   }
345 }
346
347 //____________________________________________________________
348 inline void AliTRDseedV1::SetCovRef(const Double_t *cov)
349
350 // Copy some "important" covariance matrix elements
351 //  var(y)
352 // cov(y,z)  var(z)
353 //                  var(snp)
354 //                           var(tgl)
355 //                        cov(tgl, 1/pt)  var(1/pt)
356
357   memcpy(&fRefCov[0], cov, 3*sizeof(Double_t)); // yz full covariance
358   fRefCov[3] = cov[ 5];  // snp variance 
359   fRefCov[4] = cov[ 9];  // tgl variance
360   fRefCov[5] = cov[13];  // cov(tgl, 1/pt)
361   fRefCov[6] = cov[14];  // 1/pt variance
362 }
363
364
365 //____________________________________________________________
366 inline void AliTRDseedV1::SetN(Int_t n)
367 {
368   if(n<0 || n>kNclusters) return; 
369   UInt_t mask(0x3f); 
370   fN &= ~mask; 
371   fN |= (n&mask);
372 }
373
374 //____________________________________________________________
375 inline void AliTRDseedV1::SetNUsed(Int_t n)
376 {
377   if(n<0 || n>kNclusters) return; 
378   UInt_t mask(0x3f<<6); 
379   fN &= ~mask;
380   n <<= 6; fN |= (n&mask);
381 }
382
383 //____________________________________________________________
384 inline void AliTRDseedV1::SetNShared(Int_t n)
385 {
386   if(n<0 || n>kNclusters) return; 
387   UInt_t mask(0x3f<<12); 
388   fN &= ~mask;
389   n <<= 12; fN |= (n&mask);
390 }
391
392
393 #endif
394
395
396