Add a protection againts runs with missing DCS information in the OCDB
[u/mrichter/AliRoot.git] / TRD / AliTRDcluster.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$ */
17  
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  TRD cluster                                                              //
21 //                                                                           //
22 /////////////////////////////////////////////////////////////////////////////// 
23
24 #include <TMath.h>
25
26 #include "AliLog.h"
27
28 #include "AliTRDcluster.h"
29 #include "AliTRDgeometry.h"
30 #include "AliTRDCommonParam.h"
31 #include "AliTRDtrackletWord.h"
32
33 ClassImp(AliTRDcluster)
34
35 const Int_t AliTRDcluster::fgkNlut = 128;
36 Double_t *AliTRDcluster::fgLUT = 0x0;
37
38 //___________________________________________________________________________
39 AliTRDcluster::AliTRDcluster() 
40   :AliCluster() 
41   ,fPadCol(0)
42   ,fPadRow(0)
43   ,fPadTime(0)
44   ,fLocalTimeBin(0)
45   ,fNPads(0)
46   ,fClusterMasking(0)
47   ,fDetector(0)
48   ,fQ(0)
49   ,fCenter(0)
50
51   //
52   // Default constructor
53   //
54
55   for (Int_t i = 0; i < 7; i++) {
56     fSignals[i] = 0;
57   }
58   SetBit(kLUT);
59 }
60
61 //___________________________________________________________________________
62 AliTRDcluster::AliTRDcluster(Int_t det, UChar_t col, UChar_t row, UChar_t time
63                            , const Short_t *sig, UShort_t vid) 
64   :AliCluster() 
65   ,fPadCol(col)
66   ,fPadRow(row)
67   ,fPadTime(time)
68   ,fLocalTimeBin(0)
69   ,fNPads(0)
70   ,fClusterMasking(0)
71   ,fDetector(det)
72   ,fQ(0.)
73   ,fCenter(0.)
74
75   //
76   // Constructor for self constructing cluster. In this approach the information is inserted gradualy into the 
77   // cluster and all dependencies are (re)calculated inside the cluster itself.
78   //
79   // A.Bercuci <A.Bercuci@gsi.de>
80
81   memcpy(&fSignals, sig, 7*sizeof(Short_t));
82   fQ = fSignals[2]+fSignals[3]+fSignals[4];
83   SetVolumeId(vid);
84   SetBit(kLUT);
85 }
86
87 //___________________________________________________________________________
88 AliTRDcluster::AliTRDcluster(Int_t det, Float_t q
89                            , const Float_t *pos, const Float_t *sig
90                            , const Int_t *tracks, Char_t npads, Short_t * const signals
91                            , UChar_t col, UChar_t row, UChar_t time
92                            , Char_t timebin, Float_t center, UShort_t volid)
93   :AliCluster(volid,pos[0],pos[1],pos[2],sig[0],sig[1],0.0,0x0) 
94   ,fPadCol(col)
95   ,fPadRow(row)
96   ,fPadTime(time)
97   ,fLocalTimeBin(timebin)
98   ,fNPads(npads)
99   ,fClusterMasking(0)
100   ,fDetector(det)
101   ,fQ(q)
102   ,fCenter(center)
103
104   //
105   // Constructor
106   //
107
108   for (Int_t i = 0; i < 7; i++) {
109     fSignals[i] = signals[i];
110   }
111
112   if (tracks) {
113     AddTrackIndex(tracks);
114   }
115   SetBit(kLUT);
116 }
117
118 //_____________________________________________________________________________
119 AliTRDcluster::AliTRDcluster(const AliTRDtrackletWord *const tracklet, Int_t det, UShort_t volid)
120   :AliCluster(volid,tracklet->GetX(),tracklet->GetY(),tracklet->GetZ(),0,0,0)
121   ,fPadCol(0)
122   ,fPadRow(0)
123   ,fPadTime(0)
124   ,fLocalTimeBin(0)
125   ,fNPads(0)
126   ,fClusterMasking(0)
127   ,fDetector(det)
128   ,fQ(0.)
129   ,fCenter(0.)
130 {
131   //
132   // Constructor from online tracklet 
133   //
134
135   for (Int_t i = 0; i < 7; i++) {
136     fSignals[i] = 0;
137   }
138
139 }
140
141 //_____________________________________________________________________________
142 AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
143   :AliCluster(c)
144   ,fPadCol(c.fPadCol)
145   ,fPadRow(c.fPadRow)
146   ,fPadTime(c.fPadTime)
147   ,fLocalTimeBin(c.fLocalTimeBin)
148   ,fNPads(c.fNPads)
149   ,fClusterMasking(c.fClusterMasking)
150   ,fDetector(c.fDetector)
151   ,fQ(c.fQ)
152   ,fCenter(c.fCenter)
153 {
154   //
155   // Copy constructor 
156   //
157
158   SetLabel(c.GetLabel(0),0);
159   SetLabel(c.GetLabel(1),1);
160   SetLabel(c.GetLabel(2),2);
161
162   SetY(c.GetY());
163   SetZ(c.GetZ());
164   AliCluster::SetSigmaY2(c.GetSigmaY2());
165   SetSigmaZ2(c.GetSigmaZ2());  
166
167   for (Int_t i = 0; i < 7; i++) {
168     fSignals[i] = c.fSignals[i];
169   }
170
171 }
172
173 //_____________________________________________________________________________
174 AliTRDcluster &AliTRDcluster::operator=(const AliTRDcluster &c)
175 {
176   //
177   // Assignment operator
178   //
179
180   if (&c == this) {
181      return *this;
182   }
183
184   // Call the assignment operator of the base class
185   AliCluster::operator=(c);
186
187   fPadCol         = c.fPadCol;
188   fPadRow         = c.fPadRow;
189   fPadTime        = c.fPadTime;
190   fLocalTimeBin   = c.fLocalTimeBin;
191   fNPads          = c.fNPads;
192   fClusterMasking = c.fClusterMasking;
193   fDetector       = c.fDetector;
194   fQ              = c.fQ;
195   fCenter         = c.fCenter;
196
197   SetLabel(c.GetLabel(0),0);
198   SetLabel(c.GetLabel(1),1);
199   SetLabel(c.GetLabel(2),2);
200
201   SetY(c.GetY());
202   SetZ(c.GetZ());
203   SetSigmaZ2(c.GetSigmaZ2());
204
205   for (Int_t i = 0; i < 7; i++) {
206     fSignals[i] = c.fSignals[i];
207   }
208
209   return *this;
210
211
212
213 //_____________________________________________________________________________
214 void AliTRDcluster::AddTrackIndex(const Int_t * const track)
215 {
216   //
217   // Adds track index. Currently assumed that track is an array of
218   // size 9, and up to 3 track indexes are stored in fTracks[3].
219   // Indexes are sorted according to:
220   //  1) index of max number of appearances is stored first
221   //  2) if two or more indexes appear equal number of times, the lowest
222   //     ones are stored first;
223   //
224
225   const Int_t kSize = 9;
226   Int_t  entries[kSize][2];
227
228   Int_t  i = 0;
229   Int_t  j = 0;
230   Int_t  k = 0;
231   Int_t  index;
232   Bool_t indexAdded;
233
234   for (i = 0; i < kSize; i++) {
235     entries[i][0] = -1;
236     entries[i][1] =  0;
237   }                                 
238
239   for (k = 0; k < kSize; k++) {
240
241     index      = track[k];
242     indexAdded = kFALSE; 
243
244     j = 0;
245     if (index >= 0) {
246       while ((!indexAdded) && (j < kSize)) {
247         if ((entries[j][0] == index) || 
248             (entries[j][1] ==     0)) {
249           entries[j][0] = index;
250           entries[j][1] = entries[j][1] + 1;
251           indexAdded    = kTRUE;
252         }
253         j++;
254       }
255     }
256
257   }
258
259   // Sort by number of appearances and index value
260   Int_t swap = 1;
261   Int_t tmp0;
262   Int_t tmp1;
263   while (swap > 0) {
264     swap = 0;
265     for (i = 0; i < (kSize - 1); i++) {
266       if ((entries[i][0]   >= 0) && 
267           (entries[i+1][0] >= 0)) {
268         if ((entries[i][1] < entries[i+1][1]) ||
269             ((entries[i][1] == entries[i+1][1]) &&
270              (entries[i][0] >  entries[i+1][0]))) {
271           tmp0            = entries[i][0];
272           tmp1            = entries[i][1];
273           entries[i][0]   = entries[i+1][0];
274           entries[i][1]   = entries[i+1][1];
275           entries[i+1][0] = tmp0;
276           entries[i+1][1] = tmp1;
277           swap++;
278         }
279       }
280     }
281   }               
282
283   // Set track indexes
284   for (i = 0; i < 3; i++) {
285     SetLabel(entries[i][0],i);
286   }
287
288   return;
289
290 }          
291
292 //_____________________________________________________________________________
293 void AliTRDcluster::Clear(Option_t *)
294 {
295   //
296   // Reset all member to the default value
297   //
298   fPadCol=0;
299   fPadRow=0;
300   fPadTime=0;
301   fLocalTimeBin=0;
302   fNPads=0;
303   fClusterMasking=0;
304   fDetector=0;
305   for (Int_t i=0; i < 7; i++) fSignals[i]=0;
306   fQ = 0;
307   fCenter = 0;
308   for (Int_t i = 0; i < 3; i++) SetLabel(0,i);
309   SetX(0.);
310   SetY(0.);
311   SetZ(0.);
312   AliCluster::SetSigmaY2(0.);
313   SetSigmaZ2(0.);
314   SetVolumeId(0);
315 }
316
317 //_____________________________________________________________________________
318 Float_t AliTRDcluster::GetSumS() const
319 {
320   //
321   // Returns the total charge from a not unfolded cluster
322   //
323
324   Float_t sum = 0.0;
325   for (Int_t i = 0; i < 7; i++) {
326     sum += fSignals[i];
327   }
328
329   return sum;
330
331 }
332
333 //___________________________________________________________________________
334 Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
335 {
336   //
337   // Returns the error parameterization in the radial direction for TRD clusters as function of 
338   // the calibrated time bin (tb) and optionally distance to anode wire (z). By default (no z information) 
339   // the mean value over all cluster to wire distance is chosen.
340   // 
341   // There are several contributions which are entering in the definition of the radial errors of the clusters. 
342   // Although an analytic defition should be possible for the moment this is not yet available but instead a 
343   // numerical parameterization is provided (see AliTRDclusterResolution::ProcessSigma() for the calibration 
344   // method). The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. 
345   //
346   //Begin_Html
347   //<img src="TRD/clusterXerrorDiff2D.gif">
348   //End_Html
349   //
350   // Here is a list of uncertainty components:
351   // - Time Response Function (TRF) - the major contribution. since TRF is also not symmetric (even if tail is 
352   //   cancelled) it also creates a systematic shift dependent on the charge distribution before and after the cluster.
353   // - longitudinal diffusion - increase the width of TRF and scales with square root of drift length
354   // - variation in the drift velocity within the drift cell 
355   //
356   // Author
357   // A.Bercuci <A.Bercuci@gsi.de>
358   //
359
360   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
361   static const Double_t sx[24][10]={
362     {0.000e+00, 9.352e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 2.309e+00},
363     {8.387e-02, 8.718e-02, 8.816e-02, 9.444e-02, 9.993e-02, 1.083e-01, 1.161e-01, 1.280e-01, 1.417e-01, 1.406e-01},
364     {1.097e-01, 1.105e-01, 1.127e-01, 1.151e-01, 1.186e-01, 1.223e-01, 1.272e-01, 1.323e-01, 1.389e-01, 1.490e-01},
365     {1.407e-01, 1.404e-01, 1.414e-01, 1.430e-01, 1.429e-01, 1.449e-01, 1.476e-01, 1.494e-01, 1.515e-01, 1.589e-01},
366     {1.681e-01, 1.679e-01, 1.666e-01, 1.657e-01, 1.656e-01, 1.649e-01, 1.652e-01, 1.662e-01, 1.671e-01, 1.694e-01},
367     {1.745e-01, 1.737e-01, 1.707e-01, 1.690e-01, 1.643e-01, 1.610e-01, 1.612e-01, 1.628e-01, 1.638e-01, 1.659e-01},
368     {1.583e-01, 1.558e-01, 1.535e-01, 1.488e-01, 1.445e-01, 1.419e-01, 1.428e-01, 1.451e-01, 1.462e-01, 1.494e-01},
369     {1.414e-01, 1.391e-01, 1.368e-01, 1.300e-01, 1.256e-01, 1.259e-01, 1.285e-01, 1.326e-01, 1.358e-01, 1.406e-01},
370     {1.307e-01, 1.289e-01, 1.261e-01, 1.216e-01, 1.193e-01, 1.165e-01, 1.201e-01, 1.241e-01, 1.274e-01, 1.344e-01},
371     {1.251e-01, 1.227e-01, 1.208e-01, 1.155e-01, 1.110e-01, 1.116e-01, 1.133e-01, 1.187e-01, 1.229e-01, 1.308e-01},
372     {1.234e-01, 1.209e-01, 1.175e-01, 1.127e-01, 1.094e-01, 1.093e-01, 1.109e-01, 1.155e-01, 1.210e-01, 1.275e-01},
373     {1.215e-01, 1.187e-01, 1.156e-01, 1.108e-01, 1.070e-01, 1.065e-01, 1.090e-01, 1.134e-01, 1.196e-01, 1.251e-01},
374     {1.202e-01, 1.180e-01, 1.151e-01, 1.108e-01, 1.070e-01, 1.058e-01, 1.089e-01, 1.127e-01, 1.183e-01, 1.256e-01},
375     {1.207e-01, 1.176e-01, 1.142e-01, 1.109e-01, 1.072e-01, 1.069e-01, 1.088e-01, 1.122e-01, 1.182e-01, 1.252e-01},
376     {1.213e-01, 1.182e-01, 1.156e-01, 1.102e-01, 1.076e-01, 1.063e-01, 1.091e-01, 1.132e-01, 1.181e-01, 1.243e-01},
377     {1.205e-01, 1.180e-01, 1.150e-01, 1.104e-01, 1.072e-01, 1.063e-01, 1.083e-01, 1.132e-01, 1.183e-01, 1.243e-01},
378     {1.212e-01, 1.195e-01, 1.135e-01, 1.107e-01, 1.070e-01, 1.065e-01, 1.097e-01, 1.126e-01, 1.185e-01, 1.238e-01},
379     {1.201e-01, 1.184e-01, 1.155e-01, 1.111e-01, 1.088e-01, 1.075e-01, 1.089e-01, 1.131e-01, 1.189e-01, 1.237e-01},
380     {1.197e-01, 1.186e-01, 1.147e-01, 1.113e-01, 1.085e-01, 1.077e-01, 1.105e-01, 1.137e-01, 1.188e-01, 1.245e-01},
381     {1.213e-01, 1.194e-01, 1.154e-01, 1.114e-01, 1.091e-01, 1.082e-01, 1.098e-01, 1.140e-01, 1.194e-01, 1.247e-01},
382     {1.210e-01, 1.189e-01, 1.155e-01, 1.119e-01, 1.088e-01, 1.080e-01, 1.105e-01, 1.141e-01, 1.195e-01, 1.244e-01},
383     {1.196e-01, 1.189e-01, 1.145e-01, 1.105e-01, 1.095e-01, 1.083e-01, 1.087e-01, 1.121e-01, 1.173e-01, 1.208e-01},
384     {1.123e-01, 1.129e-01, 1.108e-01, 1.110e-01, 1.080e-01, 1.065e-01, 1.056e-01, 1.066e-01, 1.071e-01, 1.095e-01},
385     {1.136e-01, 1.135e-01, 1.130e-01, 1.122e-01, 1.113e-01, 1.071e-01, 1.041e-01, 1.025e-01, 1.014e-01, 9.973e-02}
386   };
387   if(z>=0. && z<.25) return sx[tb][Int_t(z/.025)];
388   
389   Double_t m = 0.; for(Int_t id=10; id--;) m+=sx[tb][id];
390   return m*.1;
391
392 }
393
394 //___________________________________________________________________________
395 Double_t AliTRDcluster::GetSYdrift(Int_t tb, Int_t ly, Double_t/* z*/)
396 {
397   //
398   // Returns the error parameterization for TRD clusters as function of the drift length (here calibrated time bin tb)
399   // and optionally distance to anode wire (z) for the LUT r-phi cluster shape. By default (no z information) the largest 
400   // value over all cluster to wire values is chosen.
401   // 
402   // For the LUT method the dependence of s_y with x and d is obtained via a fit to the cluster to MC 
403   // resolution. (see class AliTRDclusterResolution for more details). A normalization to the reference radial position
404   // x0 = 0.675 (tb=5 for ideal vd) is also applied (see GetSYprf()). The function is *NOT* calibration aware ! 
405   // The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. A comparison 
406   // with the GAUS parameterization is also given
407   //
408   // For the GAUS method the dependence of s_y with x is *analytic* and it is expressed by the relation.
409   // BEGIN_LATEX
410   // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
411   // END_LATEX
412   // The result is displayed in the figure below as function of the drift time and compared with the LUT parameterization.
413   //Begin_Html
414   //<img src="TRD/clusterYerrorDiff2D.gif">
415   //<img src="TRD/clusterYerrorDiff1D.gif">
416   //End_Html
417   //
418   // Author
419   // A.Bercuci <A.Bercuci@gsi.de>
420   //
421
422   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
423   static const Float_t lSy[6][24] = {
424     {75.7561, 0.0325, 0.0175, 0.0174, 0.0206, 0.0232,
425      0.0253, 0.0262, 0.0265, 0.0264, 0.0266, 0.0257,
426      0.0258, 0.0261, 0.0259, 0.0253, 0.0257, 0.0261,
427      0.0255, 0.0250, 0.0259, 0.0266, 0.0278, 0.0319
428     },
429     {49.2252, 0.0371, 0.0204, 0.0189, 0.0230, 0.0261,
430      0.0281, 0.0290, 0.0292, 0.0286, 0.0277, 0.0279,
431      0.0285, 0.0281, 0.0291, 0.0281, 0.0281, 0.0282,
432      0.0272, 0.0282, 0.0282, 0.0284, 0.0310, 0.0334
433     },
434     {55.1674, 0.0388, 0.0212, 0.0200, 0.0239, 0.0271,
435      0.0288, 0.0299, 0.0306, 0.0300, 0.0296, 0.0303,
436      0.0293, 0.0290, 0.0291, 0.0294, 0.0295, 0.0290,
437      0.0293, 0.0292, 0.0292, 0.0293, 0.0316, 0.0358
438     },
439     {45.1004, 0.0411, 0.0225, 0.0215, 0.0249, 0.0281,
440      0.0301, 0.0315, 0.0320, 0.0308, 0.0318, 0.0321,
441      0.0312, 0.0311, 0.0316, 0.0315, 0.0310, 0.0308,
442      0.0313, 0.0303, 0.0314, 0.0314, 0.0324, 0.0369
443     },
444     {43.8614, 0.0420, 0.0239, 0.0224, 0.0268, 0.0296,
445      0.0322, 0.0336, 0.0333, 0.0326, 0.0321, 0.0325,
446      0.0329, 0.0326, 0.0323, 0.0322, 0.0326, 0.0320,
447      0.0329, 0.0319, 0.0314, 0.0329, 0.0341, 0.0373
448     },
449     {40.5440, 0.0434, 0.0246, 0.0236, 0.0275, 0.0311,
450      0.0332, 0.0345, 0.0347, 0.0347, 0.0340, 0.0336,
451      0.0339, 0.0344, 0.0339, 0.0341, 0.0341, 0.0342,
452      0.0345, 0.0328, 0.0341, 0.0332, 0.0356, 0.0398
453     },
454   };
455   // adjusted ...
456   return TMath::Max(lSy[ly][tb]-0.0150, 0.0010);
457
458 /*  const Double_t sy[24][10]={
459     {0.000e+00, 2.610e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 4.680e-01},
460     {3.019e-02, 3.036e-02, 3.131e-02, 3.203e-02, 3.294e-02, 3.407e-02, 3.555e-02, 3.682e-02, 3.766e-02, 3.824e-02},
461     {1.773e-02, 1.778e-02, 1.772e-02, 1.790e-02, 1.807e-02, 1.833e-02, 1.873e-02, 1.905e-02, 1.958e-02, 2.029e-02},
462     {1.774e-02, 1.772e-02, 1.746e-02, 1.738e-02, 1.756e-02, 1.756e-02, 1.739e-02, 1.720e-02, 1.743e-02, 1.769e-02},
463     {2.064e-02, 2.078e-02, 2.069e-02, 2.060e-02, 2.033e-02, 2.024e-02, 2.022e-02, 1.961e-02, 1.922e-02, 1.901e-02},
464     {2.382e-02, 2.379e-02, 2.371e-02, 2.333e-02, 2.318e-02, 2.285e-02, 2.255e-02, 2.244e-02, 2.174e-02, 2.132e-02},
465     {2.615e-02, 2.589e-02, 2.539e-02, 2.493e-02, 2.420e-02, 2.396e-02, 2.362e-02, 2.342e-02, 2.321e-02, 2.330e-02},
466     {2.640e-02, 2.638e-02, 2.577e-02, 2.548e-02, 2.477e-02, 2.436e-02, 2.416e-02, 2.401e-02, 2.399e-02, 2.402e-02},
467     {2.647e-02, 2.632e-02, 2.587e-02, 2.546e-02, 2.465e-02, 2.447e-02, 2.429e-02, 2.415e-02, 2.429e-02, 2.475e-02},
468     {2.657e-02, 2.637e-02, 2.580e-02, 2.525e-02, 2.492e-02, 2.441e-02, 2.446e-02, 2.441e-02, 2.478e-02, 2.491e-02},
469     {2.640e-02, 2.608e-02, 2.583e-02, 2.539e-02, 2.478e-02, 2.440e-02, 2.456e-02, 2.464e-02, 2.486e-02, 2.533e-02},
470     {2.636e-02, 2.630e-02, 2.584e-02, 2.542e-02, 2.483e-02, 2.451e-02, 2.449e-02, 2.467e-02, 2.496e-02, 2.554e-02},
471     {2.634e-02, 2.629e-02, 2.583e-02, 2.526e-02, 2.480e-02, 2.460e-02, 2.458e-02, 2.472e-02, 2.518e-02, 2.549e-02},
472     {2.629e-02, 2.621e-02, 2.581e-02, 2.527e-02, 2.480e-02, 2.458e-02, 2.451e-02, 2.485e-02, 2.516e-02, 2.547e-02},
473     {2.629e-02, 2.607e-02, 2.573e-02, 2.543e-02, 2.485e-02, 2.464e-02, 2.452e-02, 2.476e-02, 2.505e-02, 2.550e-02},
474     {2.635e-02, 2.613e-02, 2.578e-02, 2.523e-02, 2.491e-02, 2.465e-02, 2.470e-02, 2.467e-02, 2.515e-02, 2.564e-02},
475     {2.613e-02, 2.602e-02, 2.587e-02, 2.526e-02, 2.507e-02, 2.482e-02, 2.456e-02, 2.486e-02, 2.509e-02, 2.572e-02},
476     {2.620e-02, 2.599e-02, 2.563e-02, 2.528e-02, 2.484e-02, 2.462e-02, 2.464e-02, 2.476e-02, 2.513e-02, 2.571e-02},
477     {2.634e-02, 2.596e-02, 2.565e-02, 2.519e-02, 2.497e-02, 2.457e-02, 2.450e-02, 2.481e-02, 2.511e-02, 2.540e-02},
478     {2.593e-02, 2.589e-02, 2.563e-02, 2.511e-02, 2.472e-02, 2.453e-02, 2.452e-02, 2.474e-02, 2.501e-02, 2.543e-02},
479     {2.576e-02, 2.582e-02, 2.526e-02, 2.505e-02, 2.462e-02, 2.446e-02, 2.445e-02, 2.466e-02, 2.486e-02, 2.510e-02},
480     {2.571e-02, 2.549e-02, 2.533e-02, 2.501e-02, 2.453e-02, 2.443e-02, 2.445e-02, 2.450e-02, 2.448e-02, 2.469e-02},
481     {2.812e-02, 2.786e-02, 2.776e-02, 2.723e-02, 2.695e-02, 2.650e-02, 2.642e-02, 2.617e-02, 2.612e-02, 2.610e-02},
482     {3.251e-02, 3.267e-02, 3.223e-02, 3.183e-02, 3.125e-02, 3.106e-02, 3.067e-02, 3.010e-02, 2.936e-02, 2.927e-02}
483   };
484   if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)] -  sy[5][Int_t(z/.025)];
485
486   Double_t m = -1.e8; for(Int_t id=10; id--;) if((sy[tb][id] - sy[5][id])>m) m=sy[tb][id]-sy[5][id];
487
488   return m;*/
489 }
490
491 //___________________________________________________________________________
492 Double_t AliTRDcluster::GetSYcharge(Float_t q)
493 {
494   //
495   // Parameterization of the r-phi resolution component due to cluster charge.
496   // The value is the offset from the nominal cluster resolution defined as the
497   // cluster resolution at average cluster charge (q0).
498   // 
499   // BEGIN_LATEX
500   // #Delta #sigma_{y}(q) = a*(#frac{1}{q} - #frac{1}{q_{0}})
501   // q_{0} #approx 50
502   // END_LATEX
503   // The definition is *NOT* robust against gain fluctuations and thus two approaches are possible
504   // when residual miscalibration are available:
505   //   - determine parameterization with a resolution matching those of the gain
506   //   - define an analytic model which scales with the gain.
507   // 
508   // For more details please see AliTRDclusterResolution::ProcessCharge()
509   //
510   //Begin_Html
511   //<img src="TRD/clusterQerror.gif">
512   //End_Html
513   //
514   // Author
515   // A.Bercuci <A.Bercuci@gsi.de>
516   //
517
518   const Float_t sq0inv = 0.019962; // [1/q0]
519   const Float_t sqb    = 0.037328; // [cm]
520
521   return sqb*(1./q - sq0inv);
522 }
523
524 //___________________________________________________________________________
525 Double_t AliTRDcluster::GetSYprf(Int_t ly, Double_t center, Double_t s2)
526 {
527   //
528   // Parameterization of the cluster error in the r-phi direction due to charge sharing between 
529   // adiacent pads. Should be identical to what is provided in the OCDB as PRF [TODO]
530   //
531   // The parameterization is obtained from fitting cluster resolution at phi=exb and |x-0.675|<0.225. 
532   // For more details see AliTRDclusterResolution::ProcessCenter().
533   //
534   //Begin_Html
535   //<img src="TRD/clusterPRFerror.gif">
536   //End_Html
537   //
538   // Author
539   // A.Bercuci <A.Bercuci@gsi.de>
540   //
541
542 /*  const Float_t scy[AliTRDgeometry::kNlayer][4] = {
543     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
544     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
545     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
546     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
547     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
548     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
549   };*/
550   const Float_t lPRF[] = {0.438, 0.403, 0.393, 0.382, 0.376, 0.345};
551
552   return s2*TMath::Gaus(center, 0., lPRF[ly]);
553 }
554
555
556 //___________________________________________________________________________
557 Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
558 {
559   //
560   // Drift length correction [cm]. Due to variation of mean drift velocity along the drift region
561   // from nominal vd at xd->infinity. For drift velocity determination based on tracking information 
562   // the correction should be negligible.
563   //Begin_Html
564   //<img src="TRD/clusterXcorr.gif">
565   //End_Html
566   // TODO to be parametrized in term of drift velocity at infinite drift length
567   // A.Bercuci (Mar 28 2009)
568   //
569
570   if(tb<0 || tb>=24) return 0.;
571   const Int_t nd = 5;
572   static const Double_t dx[24][nd]={
573     {+1.747e-01,+3.195e-01,+1.641e-01,+1.607e-01,+6.002e-01},
574     {+5.468e-02,+5.760e-02,+6.365e-02,+8.003e-02,+1.067e-01},
575     {-6.327e-02,-6.339e-02,-6.423e-02,-6.900e-02,-7.949e-02},
576     {-1.417e-01,-1.424e-01,-1.450e-01,-1.465e-01,-1.514e-01},
577     {-1.637e-01,-1.619e-01,-1.622e-01,-1.613e-01,-1.648e-01},
578     {-1.386e-01,-1.334e-01,-1.261e-01,-1.276e-01,-1.314e-01},
579     {-8.799e-02,-8.299e-02,-7.861e-02,-8.038e-02,-8.436e-02},
580     {-5.139e-02,-4.849e-02,-4.641e-02,-4.965e-02,-5.286e-02},
581     {-2.927e-02,-2.773e-02,-2.807e-02,-3.021e-02,-3.378e-02},
582     {-1.380e-02,-1.229e-02,-1.335e-02,-1.547e-02,-1.984e-02},
583     {-4.168e-03,-4.601e-03,-5.462e-03,-8.164e-03,-1.035e-02},
584     {+2.044e-03,+1.889e-03,+9.603e-04,-1.342e-03,-3.736e-03},
585     {+3.568e-03,+3.581e-03,+2.391e-03,+2.942e-05,-1.585e-03},
586     {+4.403e-03,+4.571e-03,+3.509e-03,+8.703e-04,-1.425e-03},
587     {+4.941e-03,+4.808e-03,+3.284e-03,+1.105e-03,-1.208e-03},
588     {+5.124e-03,+5.022e-03,+4.305e-03,+2.023e-03,-1.145e-03},
589     {+4.882e-03,+4.008e-03,+3.408e-03,+7.886e-04,-1.356e-03},
590     {+3.852e-03,+3.539e-03,+2.057e-03,+1.670e-04,-1.993e-03},
591     {+2.154e-03,+2.111e-03,+5.723e-04,-1.254e-03,-3.256e-03},
592     {+1.755e-03,+2.101e-03,+9.516e-04,-1.649e-03,-3.394e-03},
593     {+1.617e-03,+1.662e-03,+4.169e-04,-9.843e-04,-4.309e-03},
594     {-9.204e-03,-9.069e-03,-1.182e-02,-1.458e-02,-1.880e-02},
595     {-6.727e-02,-6.820e-02,-6.804e-02,-7.134e-02,-7.615e-02},
596     {-1.802e-01,-1.733e-01,-1.633e-01,-1.601e-01,-1.632e-01}
597   };
598 //   const Double_t dx[24][nd]={
599 //     {+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00},
600 //     {-2.763e-04,-2.380e-04,-6.286e-04,-9.424e-04,+1.046e-03,+1.932e-03,+1.620e-03,+1.951e-03,-1.321e-03,-1.115e-03},
601 //     {-1.825e-03,-9.245e-04,-1.012e-03,-8.215e-04,+2.703e-05,+1.403e-03,+2.340e-03,+2.577e-03,+2.017e-03,+8.006e-04},
602 //     {-3.070e-03,-8.563e-04,-1.257e-03,+8.491e-05,+4.503e-04,-2.467e-05,-1.793e-04,+5.085e-04,+1.321e-03,+4.056e-04},
603 //     {-3.637e-03,-2.857e-03,-3.098e-03,-2.304e-03,-1.467e-03,-1.755e-03,+4.585e-04,+2.757e-03,+3.184e-03,+3.525e-03},
604 //     {-9.884e-03,-7.695e-03,-7.290e-03,-3.990e-03,-9.982e-04,+2.226e-03,+3.375e-03,+6.419e-03,+7.209e-03,+6.891e-03},
605 //     {-6.844e-03,-5.807e-03,-4.012e-03,-1.566e-03,+5.035e-04,+2.024e-03,+3.225e-03,+3.918e-03,+5.942e-03,+6.024e-03},
606 //     {-2.628e-03,-2.201e-03,-4.562e-04,+9.832e-04,+3.411e-03,+2.062e-03,+1.526e-03,+9.350e-04,+8.842e-04,+1.007e-03},
607 //     {+6.603e-04,+1.545e-03,+1.681e-03,+1.918e-03,+2.165e-03,+1.825e-03,+1.691e-03,-1.923e-04,+1.835e-04,-1.284e-03},
608 //     {+1.895e-03,+1.586e-03,+2.000e-03,+3.537e-03,+2.526e-03,+1.316e-03,+8.229e-04,-7.671e-05,-2.175e-03,-3.528e-03},
609 //     {+2.927e-03,+3.369e-03,+3.603e-03,+2.675e-03,+2.737e-03,+1.133e-03,+4.318e-04,-1.215e-03,-2.443e-03,-3.116e-03},
610 //     {+3.072e-03,+3.564e-03,+3.612e-03,+3.149e-03,+2.768e-03,+1.186e-03,+3.083e-04,-1.447e-03,-2.480e-03,-3.263e-03},
611 //     {+2.697e-03,+3.565e-03,+3.759e-03,+2.855e-03,+2.909e-03,+6.564e-04,-5.224e-04,-3.309e-04,-1.636e-03,-3.739e-03},
612 //     {+3.633e-03,+3.232e-03,+3.727e-03,+3.024e-03,+3.365e-03,+1.598e-03,-6.903e-04,-1.039e-03,-3.176e-03,-4.472e-03},
613 //     {+2.999e-03,+3.942e-03,+3.322e-03,+3.162e-03,+1.978e-03,+1.657e-03,-4.760e-04,-8.343e-04,-2.346e-03,-3.281e-03},
614 //     {+3.734e-03,+3.098e-03,+3.435e-03,+2.512e-03,+2.651e-03,+1.745e-03,+9.424e-04,-1.404e-03,-3.177e-03,-4.444e-03},
615 //     {+3.204e-03,+4.003e-03,+3.068e-03,+2.697e-03,+3.187e-03,+3.878e-04,-1.124e-04,-1.855e-03,-2.584e-03,-3.807e-03},
616 //     {+2.653e-03,+3.631e-03,+2.327e-03,+3.460e-03,+1.810e-03,+1.244e-03,-3.651e-04,-2.664e-04,-2.307e-03,-3.642e-03},
617 //     {+2.538e-03,+3.208e-03,+2.390e-03,+3.519e-03,+1.763e-03,+1.330e-04,+1.669e-04,-1.422e-03,-1.685e-03,-3.519e-03},
618 //     {+2.605e-03,+2.465e-03,+2.771e-03,+2.966e-03,+2.361e-03,+6.029e-04,-4.435e-04,-1.876e-03,-1.694e-03,-3.757e-03},
619 //     {+2.866e-03,+3.315e-03,+3.146e-03,+2.117e-03,+1.933e-03,+9.339e-04,+9.556e-04,-1.314e-03,-3.615e-03,-3.558e-03},
620 //     {+4.002e-03,+3.543e-03,+3.631e-03,+4.127e-03,+1.919e-03,-2.852e-04,-9.484e-04,-2.060e-03,-4.477e-03,-5.491e-03},
621 //     {+6.029e-03,+5.147e-03,+4.286e-03,+2.215e-03,+9.240e-04,-1.554e-03,-2.366e-03,-3.635e-03,-5.372e-03,-6.467e-03},
622 //     {+3.941e-03,+3.995e-03,+5.638e-04,-3.332e-04,-2.539e-03,-3.764e-03,-3.647e-03,-4.900e-03,-5.414e-03,-5.202e-03}
623 //   };
624   if(z>=0. && z<.25) return dx[tb][Int_t(z/.025)];
625
626   Double_t m = 0.; for(Int_t id=nd; id--;) m+=dx[tb][id];
627   return m/nd;
628 }
629
630 //___________________________________________________________________________
631 Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
632 {
633   //
634   // PRF correction for the LUT r-phi cluster shape.
635   //Begin_Html
636   //<img src="TRD/clusterYcorr.gif">
637   //End_Html
638   //
639
640   static const Float_t cy[AliTRDgeometry::kNlayer][3] = {
641     { 4.014e-04, 8.605e-03, -6.880e+00},
642     {-3.061e-04, 9.663e-03, -6.789e+00},
643     { 1.124e-03, 1.105e-02, -6.825e+00},
644     {-1.527e-03, 1.231e-02, -6.777e+00},
645     { 2.150e-03, 1.387e-02, -6.783e+00},
646     {-1.296e-03, 1.486e-02, -6.825e+00}
647   }; 
648
649   return cy[ly][0] + cy[ly][1] * TMath::Sin(cy[ly][2] * y);
650 }
651
652 //_____________________________________________________________________________
653 Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd
654                              , const Double_t *const /*q*/
655                              , const Double_t *const /*xq*/
656                              , Double_t /*z*/)
657 {
658   //
659   // (Re)Calculate cluster position in the x direction in local chamber coordinates (with respect to the anode wire 
660   // position) using all available information from tracking.
661   // Input parameters:
662   //   t0 - calibration aware trigger delay [us]
663   //   vd - drift velocity in the region of the cluster [cm/us]
664   //   z  - distance to the anode wire [cm]. By default average over the drift cell width.
665   //   q & xq - array of charges and cluster positions from previous clusters in the tracklet [a.u.]
666   // Output values :
667   //   return x position of the cluster with respect to the 
668   //   anode wire using all tracking information
669   //
670   // The estimation of the radial position is based on calculating the drift time and the drift velocity at the point of 
671   // estimation. The drift time can be estimated according to the expression:
672   // BEGIN_LATEX
673   // t_{drift} = t_{bin} - t_{0} - t_{cause}(x) - t_{TC}(q_{i-1}, q_{i-2}, ...)
674   // END_LATEX
675   // where t_0 is the delay of the trigger signal. t_cause is the causality delay between ionisation electrons hitting 
676   // the anode and the registration of maximum signal by the electronics - it is due to the rising time of the TRF 
677   // A second order correction here comes from the fact that the time spreading of charge at anode is the convolution of
678   // TRF with the diffusion and thus cross-talk between clusters before and after local clusters changes with drift length. 
679   // t_TC is the residual charge from previous (in time) clusters due to residual tails after tail cancellation. 
680   // This tends to push cluster forward and depends on the magnitude of their charge.
681   //
682   // The drift velocity varies with the drift length (and distance to anode wire) as described by cell structure simulation. 
683   // Thus one, in principle, can calculate iteratively the drift length from the expression:
684   // BEGIN_LATEX
685   // x = t_{drift}(x)*v_{drift}(x)
686   // END_LATEX
687   // In practice we use a numerical approach (AliTRDcluster::GetXcorr()) to correct for anisochronity obtained from MC 
688   // comparison (see AliTRDclusterResolution::ProcessSigma()). Also the calibration of 0 approximation (no x dependence)
689   // for t_cause is obtained from MC comparisons and impossible to disentangle in real life from trigger delay.
690   //
691   // Author
692   // Alex Bercuci <A.Bercuci@gsi.de>
693   //
694
695   AliTRDCommonParam *cp = AliTRDCommonParam::Instance(); 
696   Double_t fFreq = cp->GetSamplingFrequency();
697
698   //drift time corresponding to the center of the time bin
699   Double_t td = (fPadTime + .5)/fFreq; // [us] 
700   // correction for t0
701   td -= t0;
702   // time bin corrected for t0
703   // Bug in TMath::Nint().root-5.23.02
704   // TMath::Nint(3.5) = 4 and TMath::Nint(4.5) = 4
705   Double_t tmp = td*fFreq;
706   fLocalTimeBin = Char_t(TMath::Floor(tmp));
707   if(tmp-fLocalTimeBin > .5) fLocalTimeBin++;
708   if(td < .2) return 0.;
709   // TRF rising time (fitted)
710   // It should be absorbed by the t0. For the moment t0 is 0 for simulations.
711   // A.Bercuci (Mar 26 2009)
712   td -= 0.189;
713
714   // apply fitted correction 
715   Float_t x = td*vd + (HasXcorr() ? GetXcorr(fLocalTimeBin) : 0.);
716   if(x>0.&&x<.5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()) SetInChamber();
717
718   return x;
719 }
720
721 //_____________________________________________________________________________
722 Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W
723                              , Double_t *const y1, Double_t *const y2)
724 {
725   //
726   // Calculate, in tracking cooordinate system, the r-phi offset the cluster
727   // from the middle of the center pad. Three possible methods are implemented:
728   //   - Center of Gravity (COG) see AliTRDcluster::GetDYcog()
729   //   - Look-up Table (LUT) see AliTRDcluster::GetDYlut()
730   //   - Gauss shape (GAUS) see AliTRDcluster::GetDYgauss()
731   // In addition for the case of LUT method position corrections are also
732   // applied (see AliTRDcluster::GetYcorr())
733   //
734
735   if(IsRPhiMethod(kCOG)) GetDYcog();
736   else if(IsRPhiMethod(kLUT)) GetDYlut();
737   else if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
738   else return 0.;
739
740   if(y1) (*y1)*=W;
741   if(y2) (*y2)*=W;
742
743   return y0+fCenter*W+(IsRPhiMethod(kLUT)?GetYcorr(AliTRDgeometry::GetLayer(fDetector), fCenter):0.);
744 }
745
746 //___________________________________________________________________________
747 void AliTRDcluster::SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x, Float_t z, Float_t tgp)
748 {
749   //
750   // Set variance of TRD cluster in the r-phi direction for each method.
751   // Parameters :
752   //   - s2  - variance due to PRF width for the case of Gauss model. Replaced by parameterization in case of LUT.
753   //   - dt  - transversal diffusion coeficient 
754   //   - exb - tg of lorentz angle
755   //   - x   - drift length - with respect to the anode wire
756   //   - z   - offset from the anode wire
757   //   - tgp - local tangent of the track momentum azimuthal angle
758   //
759   // The ingredients from which the error is computed are:
760   //   - PRF (charge sharing on adjacent pads)  - see AliTRDcluster::GetSYprf()
761   //   - diffusion (dependence with drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSYdrift()
762   //   - charge of the cluster (complex dependence on gain and tail cancellation) - see AliTRDcluster::GetSYcharge()
763   //   - lorentz angle (dependence on the drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSX()
764   //   - track angle (superposition of charges on the anode wire) - see AliTRDseedV1::Fit()
765   //   - projection of radial(x) error on r-phi due to fixed value assumed in tracking for x - see AliTRDseedV1::Fit()
766   //
767   // The last 2 contributions to cluster error can be estimated only during tracking when the track angle 
768   // is known (tgp). For this reason the errors (and optional position) of TRD clusters are recalculated during 
769   // tracking and thus clusters attached to tracks might differ from bare clusters.
770   // 
771   // Taking into account all contributions one can write the the TRD cluster error parameterization as:
772   // BEGIN_LATEX
773   // #sigma_{y}^{2} = (#sigma_{diff}*Gauss(0, s_{ly}) + #delta_{#sigma}(q))^{2} + tg^{2}(#alpha_{L})*#sigma_{x}^{2} + tg^{2}(#phi-#alpha_{L})*#sigma_{x}^{2}+[tg(#phi-#alpha_{L})*tg(#alpha_{L})*x]^{2}/12
774   // END_LATEX
775   // From this formula one can deduce a that the simplest calibration method for PRF and diffusion contributions is 
776   // by measuring resolution at B=0T and phi=0. To disentangle further the two remaining contributions one has 
777   // to represent s2 as a function of drift length. 
778   // 
779   // In the gaussian model the diffusion contribution can be expressed as:
780   // BEGIN_LATEX
781   // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
782   // END_LATEX
783   // thus resulting the PRF contribution. For the case of the LUT model both contributions have to be determined from 
784   // the fit (see AliTRDclusterResolution::ProcessCenter() for details).
785   // 
786   // Author:
787   // A.Bercuci <A.Bercuci@gsi.de>
788   //
789
790   Float_t sigmaY2 = 0.;
791   Int_t ly = AliTRDgeometry::GetLayer(fDetector);
792   if(IsRPhiMethod(kCOG)) sigmaY2 = 4.e-4;
793   else if(IsRPhiMethod(kLUT)){ 
794     Float_t sd = GetSYdrift(fLocalTimeBin, ly, z); //printf("drift[%6.2f] ", 1.e4*sd);
795     sigmaY2 = GetSYprf(ly, fCenter, sd); //printf("PRF[%6.2f] ", 1.e4*sigmaY2);
796     // add charge contribution TODO scale with respect to s2
797     sigmaY2+= GetSYcharge(TMath::Abs(fQ)); //printf("Q[%6.2f] ", 1.e4*sigmaY2);
798     sigmaY2 = TMath::Max(sigmaY2, Float_t(0.0010)); //!! protection 
799     sigmaY2*= sigmaY2;
800   } else if(IsRPhiMethod(kGAUS)){
801     // PRF contribution
802     sigmaY2 = s2;
803     // Diffusion contribution
804     Double_t sD2 = dt/(1.+exb); sD2 *= sD2; sD2 *= x;
805     sigmaY2+= sD2; 
806     // add charge contribution TODO scale with respect to s2
807     //sigmaY2+= GetSYcharge(TMath::Abs(fQ));
808   }
809
810   // store tg^2(phi-a_L) and tg^2(a_L)
811   // limit parametrization to a maximum angle of 25 deg
812   if(TMath::Abs(tgp)>0.466) tgp = (tgp>0.)?0.466:-0.466;
813   Double_t tgg = (tgp-exb)/(1.+tgp*exb); tgg *= tgg;
814   Double_t exb2= exb*exb;
815
816   // Lorentz angle shift contribution 
817   Float_t sx = GetSX(fLocalTimeBin, z); sx*=sx;
818   sigmaY2+= exb2*sx; //printf("Al[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
819
820   // Radial contribution due to not measuring x in Kalman model 
821   sigmaY2+= tgg*sx; //printf("x[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
822
823   // Track angle contribution
824   sigmaY2+= tgg*x*x*exb2/12.; //printf("angle[%6.2f]\n", 1.e4*TMath::Sqrt(sigmaY2));
825
826   AliCluster::SetSigmaY2(sigmaY2);
827
828 }
829
830 //_____________________________________________________________________________
831 Bool_t AliTRDcluster::IsEqual(const TObject *o) const
832 {
833   //
834   // Compare relevant information of this cluster with another one
835   //
836   
837   const AliTRDcluster *inCluster = dynamic_cast<const AliTRDcluster*>(o);
838   if (!o || !inCluster) return kFALSE;
839
840   if ( AliCluster::GetX() != inCluster->GetX() ) return kFALSE;
841   if ( AliCluster::GetY() != inCluster->GetY() ) return kFALSE;
842   if ( AliCluster::GetZ() != inCluster->GetZ() ) return kFALSE;
843   if ( fQ != inCluster->fQ ) return kFALSE;
844   if ( fDetector != inCluster->fDetector ) return kFALSE;
845   if ( fPadCol != inCluster->fPadCol ) return kFALSE;
846   if ( fPadRow != inCluster->fPadRow ) return kFALSE;
847   if ( fPadTime != inCluster->fPadTime ) return kFALSE;
848   if ( fClusterMasking != inCluster->fClusterMasking ) return kFALSE;
849   if ( IsInChamber() != inCluster->IsInChamber() ) return kFALSE;
850   if ( IsShared() != inCluster->IsShared() ) return kFALSE;
851   if ( IsUsed() != inCluster->IsUsed() ) return kFALSE;
852   
853   return kTRUE;
854
855 }
856
857 //_____________________________________________________________________________
858 void AliTRDcluster::Print(Option_t *o) const
859 {
860   //
861   // Print cluster information
862   //
863
864   if(strcmp(o, "a")==0) {
865     Char_t mcInfo[100]; if(IsMCcluster()) snprintf(mcInfo, 100, "\n         MC[%5d %5d %5d]", GetLabel(0), GetLabel(1), GetLabel(2));
866     AliInfo(Form(
867     "\nDet[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%6.2f] Qr[%4d] FLAG[in(%c) use(%c) sh(%c)] Y[%s]"
868     "\n         LChC[c(%3d) r(%2d) t(%2d)] t-t0[%2d] Npad[%d] cen[%5.3f] mask[%d]"
869     "\n         QS[%s][%3d %3d %3d %3d %3d %3d %3d]"
870     "\n         S2Y[%e] S2Z[%e]"
871     "%s"
872     , fDetector, GetX(), GetY(), GetZ(), fQ, GetRawQ(),  
873     IsInChamber() ? 'y' : 'n', 
874     IsUsed() ? 'y' : 'n', 
875     IsShared() ? 'y' : 'n',
876     IsRPhiMethod(kGAUS)?"GAUS":(IsRPhiMethod(kLUT)?"LUT":"COG")
877     , fPadCol, fPadRow, fPadTime, fLocalTimeBin, fNPads, fCenter, fClusterMasking
878     , TestBit(kRawSignals)?"raw":"cal", fSignals[0], fSignals[1], fSignals[2], fSignals[3]
879     , fSignals[4], fSignals[5], fSignals[6]
880     , GetSigmaY2(), GetSigmaZ2()
881     , IsMCcluster()?mcInfo:""
882     ));
883   } else { 
884     AliInfo(Form("Det[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%6.2f] FLAG[in(%c) use(%c) sh(%c)] Y[%s]",
885     fDetector, GetX(), GetY(), GetZ(), fQ, 
886     IsInChamber() ? 'y' : 'n', 
887     IsUsed() ? 'y' : 'n', 
888     IsShared() ? 'y' : 'n',
889     IsRPhiMethod(kGAUS)?"GAUS":(IsRPhiMethod(kLUT)?"LUT":"COG")
890     ));
891   }
892 }
893
894 //_____________________________________________________________________________
895 void AliTRDcluster::SetPadMaskedPosition(UChar_t position)
896 {
897   //
898   // Store the pad corruption position code
899   // 
900   // Code: 1 = left cluster
901   //       2 = middle cluster;
902   //       4 = right cluster
903   //
904
905   for (Int_t ipos = 0; ipos < 3; ipos++) {
906     if (TESTBIT(position, ipos))
907       SETBIT(fClusterMasking, ipos);
908   }
909 }
910
911 //_____________________________________________________________________________
912 void AliTRDcluster::SetPadMaskedStatus(UChar_t status)
913 {
914   //
915   // Store the status of the corrupted pad
916   //
917   // Code: 2 = noisy
918   //       4 = Bridged Left
919   //       8 = Bridged Right
920   //      32 = Not Connected
921   //
922
923   for (Int_t ipos = 0; ipos < 5; ipos++) {
924     if(TESTBIT(status, ipos))
925       SETBIT(fClusterMasking, ipos + 3);
926   }
927
928 }
929
930 //___________________________________________________________________________
931 Float_t AliTRDcluster::GetDYcog(const Double_t *const, const Double_t *const)
932 {
933   //
934   // Get COG position
935   // Used for clusters with more than 3 pads - where LUT not applicable
936   //
937
938   Double_t sum = fSignals[1]
939                 +fSignals[2]
940                 +fSignals[3] 
941                 +fSignals[4]
942                 +fSignals[5];
943
944   // ???????????? CBL
945   // Go to 3 pad COG ????
946   // ???????????? CBL
947   fCenter = (0.0 * (-fSignals[1] + fSignals[5])
948                     + (-fSignals[2] + fSignals[4])) / sum;
949
950   return fCenter;
951 }
952
953 //___________________________________________________________________________
954 Float_t AliTRDcluster::GetDYlut(const Double_t *const, const Double_t *const)
955 {
956   //
957   // Calculates the cluster position using the lookup table.
958   // Method provided by Bogdan Vulpescu.
959   //
960
961   if(!fgLUT) FillLUT();
962
963   Double_t ampL = fSignals[2],
964            ampC = fSignals[3],
965            ampR = fSignals[4];
966   Int_t ilayer = AliTRDgeometry::GetLayer(fDetector);
967
968   Double_t x    = 0.0;
969   Double_t xmin, xmax, xwid;
970
971   Int_t    side = 0;
972   Int_t    ix;
973
974   Double_t xMin[AliTRDgeometry::kNlayer] = { 
975     0.006492, 0.006377, 0.006258, 0.006144, 0.006030, 0.005980 
976   };
977   Double_t xMax[AliTRDgeometry::kNlayer] = { 
978     0.960351, 0.965870, 0.970445, 0.974352, 0.977667, 0.996101 
979   };
980
981   if      (ampL > ampR) {
982     x    = (ampL - ampR) / ampC;
983     side = -1;
984   } 
985   else if (ampL < ampR) {
986     x    = (ampR - ampL) / ampC;
987     side = +1;
988   }
989
990   if (ampL != ampR) {
991
992     xmin = xMin[ilayer] + 0.000005;
993     xmax = xMax[ilayer] - 0.000005;
994     xwid = (xmax - xmin) / 127.0;
995
996     if      (x < xmin) fCenter = 0.0000;
997     else if (x > xmax) fCenter = side * 0.5000;
998     else {
999       ix      = (Int_t) ((x - xmin) / xwid);
1000       fCenter = side * fgLUT[ilayer*fgkNlut+ix];
1001     }
1002   } else fCenter = 0.0;
1003
1004   return fCenter;
1005 }
1006
1007 //___________________________________________________________________________
1008 Float_t AliTRDcluster::GetDYgauss(Double_t s2w, Double_t *const y1, Double_t *const y2)
1009 {
1010 //
1011 // (Re)Calculate cluster position in the y direction in local chamber coordinates using all available information from tracking.
1012 //
1013 // Input parameters:
1014 //   s2 - sigma of gaussian parameterization (see bellow for the exact parameterization)
1015 //   W  - pad width
1016 //   xd - drift length (with respect to the anode wire) [cm]
1017 //   wt - omega*tau = tg(a_L)
1018 // Output values :
1019 //   y1 and y2 - partial positions based on 2 pads clusters
1020 //   return y position of the cluster from all information
1021 //
1022 // Estimation of y coordinate is based on the gaussian approximation of the PRF. Thus one may
1023 // calculate the y position knowing the signals q_i-1, q_i and q_i+1 in the 3 adiacent pads by:
1024 // BEGIN_LATEX
1025 // y = #frac{1}{w_{1}+w_{2}}#[]{w_{1}#(){y_{0}-#frac{W}{2}+#frac{s^{2}}{W}ln#frac{q_{i}}{q_{i-1}}}+w_{2}#(){y_{0}+ #frac{W}{2}+#frac{s^{2}}{W}ln#frac{q_{i+1}}{q_{i}}}}
1026 // END_LATEX
1027 // where W is the pad width, y_0 is the position of the center pad and s^2 is given by
1028 // BEGIN_LATEX
1029 // s^{2} = s^{2}_{0} + s^{2}_{diff} (x,B) + #frac{tg^{2}(#phi-#alpha_{L})*l^{2}}{12}
1030 // END_LATEX
1031 // with s_0 being the PRF for 0 drift and track incidence phi equal to the lorentz angle a_L and the diffusion term 
1032 // being described by:
1033 // BEGIN_LATEX
1034 // s_{diff} (x,B) = #frac{D_{L}#sqrt{x}}{1+#(){#omega#tau}^{2}}
1035 // END_LATEX
1036 // with x being the drift length. The weights w_1 and w_2 are taken to be q_i-1^2 and q_i+1^2 respectively
1037 // 
1038 // Authors
1039 // Alex Bercuci <A.Bercuci@gsi.de>
1040 // Theodor Rascanu <trascanu@stud.uni-frankfurt.de>
1041 //
1042   Double_t w1 = fSignals[2]*fSignals[2];
1043   Double_t w2 = fSignals[4]*fSignals[4];
1044   Double_t w = w1+w2;
1045   if(w<1.){
1046     AliError("Missing side signals for cluster.");
1047     Print("a");
1048     return 0.;
1049   }  
1050
1051   //Double_t s2w = s2/W/W;
1052   Float_t y1r  = fSignals[2]>0 ? (-0.5 + s2w*TMath::Log(fSignals[3]/(Float_t)fSignals[2])) : 0.;
1053   Float_t y2r  = fSignals[4]>0 ? (0.5 + s2w*TMath::Log(fSignals[4]/(Float_t)fSignals[3])) : 0.;
1054
1055   if(y1) (*y1) = y1r;
1056   if(y2) (*y2) = y2r;
1057
1058   return fCenter      = (w1*y1r+w2*y2r)/w;
1059 }
1060
1061
1062
1063 //_____________________________________________________________________________
1064 void AliTRDcluster::FillLUT()
1065 {
1066   //
1067   // Create the LUT
1068   //
1069
1070   // The lookup table from Bogdan
1071   Float_t lut[AliTRDgeometry::kNlayer][fgkNlut] = {  
1072     {
1073       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611, 
1074       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187, 
1075       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704, 
1076       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160, 
1077       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557, 
1078       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901, 
1079       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207, 
1080       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483, 
1081       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727, 
1082       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952, 
1083       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157, 
1084       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345, 
1085       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520, 
1086       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683, 
1087       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824, 
1088       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978 
1089     },
1090     {
1091       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642, 
1092       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241, 
1093       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770, 
1094       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229, 
1095       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627, 
1096       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968, 
1097       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266, 
1098       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536, 
1099       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772, 
1100       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991, 
1101       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187, 
1102       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369, 
1103       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538, 
1104       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694, 
1105       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832, 
1106       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1107     },
1108     {
1109       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674, 
1110       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296, 
1111       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839, 
1112       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299, 
1113       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693, 
1114       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033, 
1115       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323, 
1116       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586, 
1117       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815, 
1118       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027, 
1119       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217, 
1120       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393, 
1121       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555, 
1122       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705, 
1123       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839, 
1124       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1125     },
1126     {
1127       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708, 
1128       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356, 
1129       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910, 
1130       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371, 
1131       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759, 
1132       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095, 
1133       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378, 
1134       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633, 
1135       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857, 
1136       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061, 
1137       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244, 
1138       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414, 
1139       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571, 
1140       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715, 
1141       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846, 
1142       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1143     },
1144     {
1145       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744, 
1146       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419, 
1147       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984, 
1148       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444, 
1149       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824, 
1150       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152, 
1151       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432, 
1152       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677, 
1153       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896, 
1154       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094, 
1155       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270, 
1156       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435, 
1157       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586, 
1158       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725, 
1159       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852, 
1160       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1161     },
1162     {
1163       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792, 
1164       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505, 
1165       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081, 
1166       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541, 
1167       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917, 
1168       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235, 
1169       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511, 
1170       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748, 
1171       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962, 
1172       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152, 
1173       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324, 
1174       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483, 
1175       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629, 
1176       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759, 
1177       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859, 
1178       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1179     }
1180   }; 
1181
1182   if(!fgLUT) fgLUT = new Double_t[AliTRDgeometry::kNlayer*fgkNlut];
1183
1184   for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1185     for (Int_t ilut  = 0; ilut  < fgkNlut; ilut++  ) {
1186       fgLUT[ilayer*fgkNlut+ilut] = lut[ilayer][ilut];
1187     }
1188   }
1189 }
1190